7

In my C++ project, I have three classes, Particle, Contact, and Network. The Network class will have N particles (std::vector<Particle> particles) and Nc contacts (std::vector<Contact> contacts). Particle objects will each have a number of contacts, represented by Contact objects, and contacts are shared by pairs of particles.

It will be necessary for me to know what contacts a specific particle has, and I have decided to be able to obtain that information through the Network object (e.g., using network.getContacts(Particle particle)).

My question can be broken into two parts:

First, what would be an efficient and community-recommended way to store the contacts for each particle?

I believe that I can either create a vector of vector of shared pointers to Contact objects, vector<vector<shared_ptr<Contact>>> prtContacts (where I have left off std:: for ease of viewing), so that prtContacts[i] will contain a vector of shared pointers to the contacts of the ith particle. It has been suggested that using shared_ptr in this context is useful, as it will ensure that I am not creating duplicate Contact objects for two particles with a shared contact. Is there an alternative method that would make more sense?

Second, if I do end up creating prtContacts as defined above, what is the best way to initialize it? My current thinking is to have an initial loop that initializes std::vector<Contact> contacts with the Nc contacts and then have a second loop that initializes prtContacts with shared pointers to the Contacts. Is that reasonable or, again, are there alternative, more efficient, approaches?

Thank you!


Some additional information:

The number of particles, N, will range from O(10^1) to perhaps O(10^4). The number of contacts per particle will most likely not exceed 6-7, but will be allowed to vary. Each contact will have a position and contact force (2-component vector) associated with it, and both positions and contact forces will be allowed to vary during the course of simulations (e.g., Monte Carlo simulations).

  • What is your typical access pattern for the most compute intensive steps of the algorithms you're applying to this data structure? – Omnifarious Nov 10 '17 at 18:53
  • @Omnifarious I added some additional information to the question that is hopefully helpful. But I will want to use this for various simulations (e.g., Monte Carlo) where I am changing certain positions and/or contact forces of the contacts or positions of the particles at every iteration. – AnInquiringMind Nov 10 '17 at 21:24
  • Will the simulation engine basically loop through the contacts and use the information to update the associated particles? Or will you be looping through the particles? The thing you loop through should be allocated in all contiguous memory. Basically as members of a vector or something like that. – Omnifarious Nov 10 '17 at 21:36
  • @Omnifarious The Network class will ideally be flexible and allow for both cases (depending on which is the main control parameter – particle positions or contact forces). That's why I have defined both as vectors. But in terms of having a list of contacts per particle, knowing the extra information I have provided, is my implementation outlined in my question okay (with shared pointers being replaced with regular pointers, as you suggested)? Or is more information still needed? Thanks for the continued help! – AnInquiringMind Nov 10 '17 at 21:44
  • 1
    @Omnifarious Then I think the best way forward might be to assign a unique contact ID/index (or a unique pair corresponding to the particle IDs) to each contact so that it doesn't matter where it is located in the vector... Though that might result in always having to rely on if statements to find the correct contact? – AnInquiringMind Nov 10 '17 at 22:45
7

The right solution will depend a lot on your non-functional requirements. How many particles will you typically have? How many contacts does a particle typically have? How much is a network updated, and in what ways? How fast does it need to be?

I suggest referring to objects by index if you can and then particle contacts can be a list of contact indexes. If a particle typically has a small number of contacts then instead of a vector I suggest using something that is optimized for a small number of elements like boost::container::small_vector, absl::inlined_vector, folly::small_vector or llvm::SmallVector.

using ParticleContacts = boost::container::small_vector<int, 12>;

class Network {
  std::vector<Particle> particles;
  std::vector<Contact> contacts;
  std::vector<ParticleContacts> particle_contacts;
 public:
  ParticleContacts getParticleContacts(int particle_index) const {
    return particle_contacts.at(particle_index);
  }
  Contact getContact(int contact_index) const {
    return contacts.at(contact_index);
  }
  //...
};

Live demo.

Edit: If you have to be able to delete from the middle of the vector of contacts then using an index will not work. I would seriously consider if you can change your algorithm so that this is not necessary. Not only does it invalidate pointers, iterators and indexes but it also is quite inefficient because all the contacts after the deletion have to be moved. Perhaps, for example, you can just leave orphaned contacts until the end of the simulation with no connection to a particle. Anyway, if you really need to be able to delete from the middle of the list of contacts perhaps a vector is not the best data structure and you should prefer something like a linked list. The advantage of a linked list is you can do quick deletion from the middle and also pointers to contacts in the list are not invalidated by inserts or deletions so your list of particle contacts can be raw pointers. The downside is it is much slower and does not allow random access.

using ParticleContacts = boost::container::small_vector<Contact*, 12>;

class Network {
  std::vector<Particle> particles;
  std::forward_list<Contact> contacts;
  std::vector<ParticleContacts> particle_contacts;
 public:
  ParticleContacts getParticleContacts(int particle_index) const {
    return particle_contacts.at(particle_index);
  } 
  void addParticle(const Particle& particle) {
    particles.push_back(particle);
    particle_contacts.emplace_back();
  } 
  Contact* addParticleContact(int particle_index_1, int particle_index_2) {
    contacts.emplace_front();
    particle_contacts.at(particle_index_1).push_back(&contacts.front());
    particle_contacts.at(particle_index_2).push_back(&contacts.front());  
    return &contacts.front();
  }
};

Live demo.

Alternatively, you could use std::vector<std::unique_ptr<Contact>>. Raw pointers to the contact won't be invalidated by insertions and deletions. It has the advantage of allowing random access but the disadvantage of being slower at insertions and deletions: Live demo

The right choice will depend on how the data is typically accessed by the simulation.

  • Thanks for the response! I updated my initial question with some more information – please let me know if you need anything else. The one thing with contact indices is that won't they change as contacts are added and removed (the number of contacts, Nc, is not necessarily fixed). In that case, aren't pointers better than indices? – AnInquiringMind Nov 10 '17 at 22:25
  • Unless you mean a unique contact index that I assign to each contact, and not the index in the vector... – AnInquiringMind Nov 10 '17 at 22:31
  • 1
    You only need to use shared_ptr if you need shared ownership. If I understand things correctly the network can have unique ownership of the contacts and then the list of particle contacts can be non owning references like raw pointers. But even unique_ptr might be too much overhead. – Chris Drew Nov 11 '17 at 7:22
  • 1
    @PhysicsCodingEnthusiast: an adjacency matrix would only make sense here when implemented as a sparse matrix (since I guess you do not want to reserve memory for 10^8 contacts if 10^4 * 7 would be enough, I guess?) But since sparse matrices are implemented typically with row or column lists, this leads effectively to the same solution shown in this answer. – Doc Brown Nov 11 '17 at 19:22
  • 1
    ... using unique IDs is another possibility here. If you go that route, use an ID for particles as well as for contacts and use a hashmap container for both, where the ID is the key. That will allow fast insert and deletion as well. – Doc Brown Nov 11 '17 at 19:29
3

So, this really depends on the access patterns you anticipate. First, shared_ptr comes at a potentially large cost in a multi-threaded program. Incrementing and decrementing the reference counts will cause processor cache invalidation on other CPUs. This will cause those operations to be slightly more expensive, and it will cause a lot of cache thrashing if those objects are being accessed in any way from multiple threads because even read access to bytes near the reference count will result in a very costly main memory fetch after the reference count has been updated.

Now, one way you could mitigate this is to not use make_shared like you're told to. This will result in the reference count being allocated separately from the object. But that potentially moves the cache problem elsewhere.

The various kinds of pointer are about expressing ownership. Does the object really have shared ownership? Or is it more appropriate to think of the Network as a whole as owning the object? If it's the latter, it might make more sense to store bare pointers and simply have a few vector objects (one per type) in the Network that hold unique_ptrs to all of the graph elements.

Is the number of contacts per particle bounded at a reasonably small number? If so, it might make sense to use a fixed size array that stores a list of bare pointers to contacts and have a Network object that actually owns all of the objects.

I'm sure there are more considerations. And my primary focus here has been on performance rather than some abstract notion of the 'right' way to do things. Perhaps something other than performance is your primary goal. Though, if that's the case, why are you using C++ here?

  • Thanks for the response! I can start responding here; if there is too much information, please let me know, and I can add the relevant portions to my question. The Network object should own all particles and contacts, and I will have a vector of Nc Contact objects called contacts (see my question). But if contacts[i] is between particles[j] and particles[k], and I want that information to be obtainable by network.getContacts(particles[j]) and network.getContacts(particles[k]), won't that require two pointers to the same object (which is the purpose of shared_ptr, right?)? – AnInquiringMind Nov 10 '17 at 19:02
  • The number of contacts per particle should never exceed 6-7, but the number is not fixed (so a vector is preferable). I definitely care about performance, which I usually take to mean "right" in the context of C++! – AnInquiringMind Nov 10 '17 at 19:07
  • No, shared_ptr should not be your goto whenever you have two pointers pointing to the same thing. It's about lifetime management. If you can't think of a sane way to manage the lifetime of an object than to have shared ownership semantics, use shared_ptr. Otherwise, don't. For example, I use a reference counted pointer in something where I have buffers that can be built out of fragments from other buffers. I don't copy any parts of the buffer, I refer to them from buffers that logically contain them. shared_ptr is the only way there. – Omnifarious Nov 10 '17 at 19:08
  • Okay, so suppose that everything in my post is the same, except that I use regular pointers rather than shared pointers. Do my approaches to both questions make sense, or are there more efficient/performant approaches? Next, if I use regular pointers, should I have two pointers to the same Contact object, or should I try to use the same pointer for the vector of contacts for both particles (if that's even possible)? Thanks! – AnInquiringMind Nov 10 '17 at 19:21
  • @PhysicsCodingEnthusiast - Two pointers to the same Contact. Leave it up to the Network to make sure that the pointers get deleted from the appropriate Particles before the Contact is deleted. – Omnifarious Nov 10 '17 at 20:12
1

I'd largely echo Chris' answer but with some possible deviations. The ultimate thing I'd echo wholeheartedly is to store indices to your particles and contacts, stored in big sequences like a couple of big vector pers network.

However, one part where I disagree is that you can use indices and a random-access sequence and still have constant-time removals (with a super small constant) without index invalidation if you use an indexed free list/holes strategy. If particles die, then you can remove them in constant-time and reclaim/insert in constant-time without using any extra memory, like so:

enter image description here

... where union Nodes nodes[n] might be std::vector<union Node>. It takes a little more elbow grease if Particle or Contact isn't trivially constructible or destructible (making a proper container in C++ for non-PODs is a bit of a PITA with placement new, aligned storage, and manual dtor invocations). Hopefully they are, as generally can be the case for this type of data.

First, what would be an efficient and community-recommended way to store the contacts for each particle?

Assuming it's efficient to just simulate and gather the contacts for all particles at once (which is often the case in many simulations), short answer:

// Stores all the contacts between particles. Each per-particle 
// contact list begins with a count followed by the indices to each
// contact.
std::vector<int> contact_data;

// Stores an index into the contact data for each particle.
std::vector<int> contacts;

To get the contacts for the nth particle, we do:

// Fetch the index into the contact data for the nth particle.
const int cd_index = contacts[n];

// Fetch the number of contacts.
const int num = contact_data[cd_index];

// Fetch a pointer to the array of contact indices.
const int* indices = contact_data.data() + cd_index + 1;

And now you have the indices of the contacts and the number of them for the nth particle. You can clear the list and compute it for each time step. You can also do it multithreaded with a bit of work (gather the contacts for ranges of particles in each thread locally, then merge the results into the two vectors at the end).

network.getContacts(Particle particle)

If you follow the suggestion to use indices, the above turns to this:

// Returns the number of contacts and a pointer to the array
// of contact indices for the particle indexed by 'particle_index'.
std::pair<int, const int*> network.getContacts(int particle_index)

... you can use something a bit nicer than std::pair, like your own structure or class.

SmallVectors and the Like

Now using something like boost::small_vector or llvm::SmallVector is a very good solution, far superior to requiring a heap allocation for every particle as would be the case if you used std::vector for each one. Chris' answer is already really good. But if you want to persistently store the contacts, these small buffer optimizations start to get a little bit explosive in memory use which can translate to extra cache misses (they are pretty great for short-lived storage though on the stack), which is why I recommend just a big old vector of integers instead with the assumption that your simulation might be better off computing all the contacts at once rather than computing them on demand in getContacts (at which point a small vector would start to become much more appropriate).

[...] in this context is useful, as it will ensure that I am not creating duplicate Contact objects for two particles with a shared contact. Is there an alternative method that would make more sense?

shared_ptr is for sharing ownership of a resource and extending its lifetime until all owners release their ownership. To just refer to a resource in multiple places without sharing ownership (which is all you need here) and without duplication of the resource's data, just use plain old indices or pointers. There's a pretty hefty relative overhead to shared_ptr used on a per-particle or per-contact basis.

A quick thing in practice using the above techniques (including the indexed free list), and not really "particles" in this case but still agents that collide with each other and a collision list between agents (gathered all at once into two vectors of integers as shown above): 500,000 agents bouncing off of each other, single-threaded, with a good portion of the time spent just plotting pixels and not merely simulating:

enter image description here

First, what would be an efficient and community-recommended way to store the contacts for each particle?

If you want to do things very efficiently when the temptation is to store a boatload of tiny lists, ideally store all the data in big sequences, not a bunch of teeny ones. You don't actually need to instantiate a million containers in order to represent the equivalent of a million containers associated to a million elements. You can just use two, for example: one storing all the data and one parallel to the million elements storing starting indices into that data.

This is a general optimization strategy regardless of context, and it addresses a common performance gotcha in languages that provide a lot of convenient containers. Those containers are very efficient for storing a million elements in one container, but not very efficient for storing a million containers with a few elements each. That applies even to the ones that use SBO (small buffer optimizations) like SmallVector or std::string, because the "small buffer" is either too large (at which point we're wasting memory and getting huge strides from one element to the next) or too small and incurring heap allocations. It's generally not "just right" as is the case if you use one big buffer for all the data.

Note that this answer might be a bit overkill. I'm not really accounting for productivity at all and assuming you have a bit of extra time to deal with the inconvenience of storing everything in one big container and possibly implementing that indexed free list above (I'm not a good person to ask for efficient solutions that can be implemented in 10 minutes... well, I think I can implement these in 10 minutes but only because I've had lots of practice). For VFX it's not uncommon to deal with hundreds of millions of particles, and so I'm often required to use these types of techniques and handroll containers at the bare minimum just to get acceptable performance that won't have users complaining that some other software is faster or having studios switch to simulating with something else on their farms. But I get kinda excited to share some very general techniques to speed things up, and a good one that doesn't get too fancy and hopefully isn't too overkill is just index your elements and store everything in big containers, not a boatload of small ones. I'm actually restraining myself a lot for this answer because there's a lot more I can cover on this topic including an alternative to that indexed free list for the largest scale inputs (far more complex to implement), but I'll stop here!

0

I have to ask - why are the Contacts separate from the Particles? It seems like each Particle should contain a list of its contacts. Or if that's not appropriate, maybe another class or struct that contains a Particle and an array (or vector or whatever) of Contacts. Something like:

struct ParticleCollision {
    Particle part;
    std::vector<Contact> contacts;
};

The Network would then contain a std::vector of ParticleCollisions.

  • That's an interesting suggestion. I wonder if there is any performance concerns with doing such a successive "encapsulation". As for Particle objects containing a list of its contacts, the point is that the Network object entirely decides the relationship between particles and contacts, so all the related data should be accessible solely through the Network object. – AnInquiringMind Nov 11 '17 at 5:31
  • We already discussed this in the OPs former question - no single Particle is owner of one Contact object, since a contact is associated to two different particles. So when requesting the set of contacts of those two particles, the intersection of those sets should contain the common Contact object. – Doc Brown Nov 11 '17 at 8:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.