Make Vector/Point more efficient, drop the coordinate system pointer

Corsika operates mostly in a single coordinate system and may generate large numbers of particles with positions. In the current design, each particle with a position carries a coordinate system pointer with that position.

If we build a particle record from position and momentum of a particle, and we implement those with Point and Vector, then we carry two coordinate system pointers around, although position and momentum will usually be defined in the same coordinate system. If we consider four-vectors (see #62 (closed)), then we also would need to carry a inertial frame pointer around.

This design is inefficient. I propose to store all points and vectors and four vectors in the same internal coordinate system. When all objects are always stored in the same coordinate system, we can skip the pointer to the coordinate system. Conversions from the internal coordinate system to an external coordinate system only happens in the getters and setters of the components.

This approach scales much better, one can easily define a Particle as a struct of momentum and position without inefficiency or worry that momentum and position are defined in different coordinate systems.

changed title from Make Vector more efficient to Make Vector/Point more efficient, drop the coordinate system pointer

See also #36 (closed) (geometry-store-references-to-coordinatesystem-in-low-level-vs-high-level-objects) which is a variant of this issue in a less general form.

added Development Feature request labels

changed milestone to %Physics demonstrator (hadron core)

I am not sure why this wheel is getting reinvented again. It is clear that using Offline geometry is inefficient. In Auger case it made sense since there were gazillions of coordinate systems defined by different detector components and the manual bookkeeping would be a nightmare.

In case of Corsika there is very clear where is what. Particles live mostly in the stack where they are in the shower coordinate system. There is clear boundary where they have to be transformed, eg when sent to the hadronic models or such, and when environment is probed. Instead of coordinate-system reference in each geometrical object we need the boundary code where the transformations are done. Neither models nor environment should know anything about the original coordinate system of the particles, they get them served already in the correct representation.

In short: coordinate-system pointer in all objects --> boundary transformation code.

I think we agree on the requirements. I am arguing that since most particles will be in the same coordinate system, we do not need store a coordinate system pointer for each point and vector. The coordinate system pointer is strictly only needed when you set or get coordinates from the abstract point or vector. As a user, you don't know or need to know what the point or vector uses internally as a coordinate system.

We can exploit this to make it always the same internal coordinate system, which we then do not have to track with a pointer for each point and vector, which is very inefficient for a code that is supposed to handle millions of points and vectors.

that's why i originally opposed that you port the offline geometry directly to corsika. corsika should use only eigen classes. one can then think of how to elegantly implement the transformations, but this is all external to the low-level objects. one can have somekind of coordinate-system reference counting (ala offline coordinate-system registry) but it's not embedded into each geometry object.

I agree with this :)

Here is the summary from the meeting today: we will remove the local CoordinateSystem reference from all low-level object. The internal representation is always in a single frame (root-CS), and transformations are performed in "set" and "get" routines using provided CoordinateSystem objects.

Users can always retrieve the components, or the internal vector (QuantityVector as wrapper around eigen3) to perform local computations. This is already possible now.

The CoordinateSytem object needs to keep an internal TransformationMatrix cached for all conversion to and from the root-CS to itself.

We can have object counting to keep CoordinateSystem objects alive in a "manager" until they are not used further.

It is up to investigation if:

• relying on a single root-CS leads to more transformation and thus additional computing time
• Numerical precision far from the origin of root-CS may become a problem. Theoretically, we can in the future allow to create a fixed number (at compile time) of root-CSs relative to each other and distinguish between them with template magic. This is not clear yet how to do, and it first needs to be investigated if this will solve any upcoming problem.

@rulrich I believe all of this has been implemented but it's been many months since I looked at the geometry in detail. Can you confirm?

No, currently still each Point, Vector or other object has a shared_ptr to a CoordinateSystem, which is a sizable resource penalty. This may cost as much of 10% of computing right now.....

I don't think to drop the idea of CoordinateSystems is the solution (it would 'fix' the problem, but it would kill the feature), so a real solution is still to be defined. @lukas proposed a static version and I think this may be the way to go, but we have not seen a demonstrator so far.

Excellent! Some technical suggestion for the CoordinateSystem manager: it could use a Boost.Pool internally: https://www.boost.org/doc/libs/1_69_0/libs/pool/doc/html/boost_pool/pool/interfaces.html#boost_pool.pool.interfaces.interfaces

A pool optimises frequent allocation and deallocation of objects of the same size. It allocates a chunk of memory from the heap and then serves allocation requests from the preallocated chunk. Using malloc and free on the Pool is very fast, about 1 or 2 orders of magnitude faster than from the heap.

The library also has a singleton pool, which may be useful in this context.

I don't think we need singleton. If we pass around the environment this is probably a better solution. In the end the pool is an implementation detail, but certainly a smart one.

added Critical label

For a final assessment, we should make a numerical evaluation of the potential impact of single-frame vs. multi-frame situation. With some tests we can find out in what situations single-frame will start to yield lower performance wrt. to multi-scale. Then we can discuss those findings.

Let me emphasize that this proposal is not only about saving memory, it also simplifies the code a lot and so reduces maintenance and compile-time costs.

For example, often you will deal with collections of points in the same coordinate system. E.g. you will want a vector

std::vector<Point>

If each point stores the CS pointer, this is inefficient and someone will write a special vector which stores the pointer only once. This class has to be implemented by hand.

class SpecialEfficientVectorOfPoints

Now the next person cannot use a vector, but needs a deque or a map, so you run into the problem again. Do you want to implement also

class SpecialEfficientDequeOfPoints;
class SpecialEfficientMapOfPoints;

And what about unordered_map. And so on. A framework is most powerful if the elements are composable and work well with the C++ stdlib. I think this is even more important than the performance.

Do we have a curved atmosphere option available for testing already? A flat earth/flat atmosphere is a simpler geometry...

I haven't had the time yet to look into the geometry in corsika...

changed milestone to %ICRC2019

changed milestone to %Full physics demonstratror

added Performance label

removed Critical label

I removed the "critical" label, since obviously it is not critical after such a long time.

However, this is still an important performance-related issue and must be addressed. We cannot keep the dynamic CoordinateSystem approach as it is, since it currently creates significant overhead.

A static version of CoordinateSystems would help to resolve this -- but the solution remains unclear right now.

marked this issue as related to #703