@@ -316,30 +316,23 @@ vector< unordered_set<int> > OverlapCalculator::neighborhoods() const
316316 typedef unordered_set<int > SiteSet;
317317 typedef std::vector< boost::shared_ptr<SiteSet> > SiteSetPointers;
318318 SiteSetPointers rvptr (cntsites);
319- // helper function that initializes neighbor set at site i
320- auto initsiteset = [&](int i, int & counter) {
321- if (!rvptr[i]) {
322- rvptr[i].reset (new SiteSet (&i, &i + 1 ));
323- ++counter;
324- }
325- };
326- // create self-neighborhoods unless prohibited by mask
327- for (int j0 = 0 , cnt = 0 ; j0 < cntsites && cnt < cntsites; ++j0)
328- {
329- for (int j1 = j0; j1 < cntsites; ++j1)
330- {
331- if (!this ->getPairMask (j0, j1)) continue ;
332- initsiteset (j0, cnt);
333- initsiteset (j1, cnt);
334- }
335- }
336319 int n = this ->count ();
337320 for (int index = 0 ; index < n; ++index)
338321 {
339322 double olp = this ->suboverlap (index);
340323 if (olp <= 0.0 ) continue ;
341324 int j0 = int (this ->subvalue (SITE0_OFFSET, index));
342325 int j1 = int (this ->subvalue (SITE1_OFFSET, index));
326+ if (!rvptr[j0].get ())
327+ {
328+ rvptr[j0].reset (new SiteSet);
329+ rvptr[j0]->insert (j0);
330+ }
331+ if (!rvptr[j1].get ())
332+ {
333+ rvptr[j0]->insert (j1);
334+ rvptr[j1] = rvptr[j0];
335+ }
343336 assert (rvptr[j0].get () && rvptr[j1].get ());
344337 if (rvptr[j0].get () == rvptr[j1].get ()) continue ;
345338 // here we need to unify the j0 and j1 sets
@@ -351,6 +344,27 @@ vector< unordered_set<int> > OverlapCalculator::neighborhoods() const
351344 rvptr[*idx1] = rvptr[j0];
352345 }
353346 }
347+ // helper function that initializes neighbor set at site i
348+ auto initsiteset = [&](int i, int & counter) {
349+ if (!rvptr[i]) {
350+ rvptr[i].reset (new SiteSet (&i, &i + 1 ));
351+ ++counter;
352+ }
353+ };
354+ // count initialized items in rvptr
355+ int cnt = count_if (rvptr.begin (), rvptr.end (),
356+ [](SiteSetPointers::value_type& p) { return bool (p); });
357+ // create self-neighborhoods unless prohibited by mask
358+ for (int j0 = 0 ; j0 < cntsites && cnt < cntsites; ++j0)
359+ {
360+ for (int j1 = j0; j1 < cntsites; ++j1)
361+ {
362+ if (rvptr[j0] && rvptr[j1]) continue ;
363+ if (!this ->getPairMask (j0, j1)) continue ;
364+ initsiteset (j0, cnt);
365+ initsiteset (j1, cnt);
366+ }
367+ }
354368 unordered_set<const SiteSet*> duplicate;
355369 vector< unordered_set<int > > rv;
356370 SiteSetPointers::const_iterator ssp;
0 commit comments