Browse code

single sized buckets working

Tom Sherman authored on 22/08/2018 17:56:50
Showing5 changed files

... ...
@@ -1,154 +1,579 @@
1
-#include "GapsAssert.h"
2 1
 #include "AtomicDomain.h"
3
-#include "math/Random.h"
2
+#include "GapsAssert.h"
4 3
 
5
-#include <stdint.h>
6
-#include <utility>
4
+#include <vector>
7 5
 
8
-// O(1)
9
-Atom AtomicDomain::front() const
6
+// should only ever use this when copying a default constructed bucket
7
+AtomBucket& AtomBucket::operator=(const AtomBucket& other)
10 8
 {
11
-    return mAtoms[mAtomPositions.begin()->second];
9
+    GAPS_ASSERT(!other.mFull);
10
+    GAPS_ASSERT(other.mPrev == NULL);
11
+    GAPS_ASSERT(other.mNext == NULL);
12
+    GAPS_ASSERT(other.mOverflow == NULL);
13
+
14
+    mFull = false;
15
+    mPrev = NULL;
16
+    mNext = NULL;
17
+    mOverflow = NULL;
12 18
 }
13 19
 
14
-// O(1)
15
-Atom AtomicDomain::randomAtom() const
20
+/////////////////////// ALLOCATOR FOR OVERFLOW BUCKETS /////////////////////////
21
+
22
+static const unsigned poolsize = 1024;
23
+
24
+class AtomBucketPool
16 25
 {
17
-    uint64_t ndx = mRng.uniform64(0, mAtoms.size() - 1);
18
-    return mAtoms[ndx];
26
+public:
27
+
28
+    AtomBucketPool()
29
+        : mPool(std::vector<AtomBucket>(poolsize, AtomBucket())), mCurrent(0)
30
+    {}
31
+
32
+    AtomBucket* create()
33
+    {
34
+        return &(mPool[mCurrent++]);
35
+    }
36
+
37
+    bool depleted()
38
+    {
39
+        return mCurrent == poolsize;
40
+    }
41
+
42
+private:
43
+
44
+    std::vector<AtomBucket> mPool;
45
+    unsigned mCurrent;
46
+};
47
+
48
+class AtomBucketAllocator
49
+{
50
+public:
51
+
52
+    AtomBucketAllocator()
53
+    {
54
+        mAllPools.push_back(new AtomBucketPool());
55
+    }
56
+
57
+    ~AtomBucketAllocator()
58
+    {
59
+        for (unsigned i = 0; i < mAllPools.size(); ++i)
60
+        {
61
+            delete mAllPools[i];
62
+        }
63
+    }
64
+
65
+    AtomBucket* create()
66
+    {
67
+        if (mAllPools.back()->depleted())
68
+        {
69
+            mAllPools.push_back(new AtomBucketPool());
70
+        }
71
+        return mAllPools.back()->create();
72
+    }
73
+
74
+private:
75
+
76
+    std::vector<AtomBucketPool*> mAllPools;
77
+};
78
+
79
+static AtomBucket* createAtomBucket()
80
+{
81
+    //static AtomBucketAllocator allocator;
82
+    //return allocator.create();
83
+    return new AtomBucket();
19 84
 }
20 85
 
21
-// Average Case O(1)
22
-uint64_t AtomicDomain::randomFreePosition() const
86
+////////////////////////////////// ATOM ////////////////////////////////////////
87
+
88
+Atom::Atom() : pos(0), mass(0.f) {}
89
+
90
+Atom::Atom(uint64_t p, float m) : pos(p), mass(m) {}
91
+
92
+void Atom::operator=(Atom other)
23 93
 {
24
-    uint64_t pos = 0;
25
-    do
94
+    pos = other.pos;
95
+    mass = other.mass;
96
+}
97
+
98
+Archive& operator<<(Archive &ar, Atom &a)
99
+{
100
+    ar << a.pos << a.mass;
101
+    return ar;
102
+}
103
+
104
+Archive& operator>>(Archive &ar, Atom &a)
105
+{
106
+    ar >> a.pos >> a.mass;
107
+    return ar;
108
+}
109
+
110
+//////////////////////////// ATOM NEIGHBORHOOD /////////////////////////////////
111
+
112
+
113
+AtomNeighborhood::AtomNeighborhood(Atom *l, Atom *c, Atom *r)
114
+    : center(c), left(l), right(r)
115
+{}
116
+
117
+bool AtomNeighborhood::hasLeft()
118
+{
119
+    return left != NULL;
120
+}
121
+
122
+bool AtomNeighborhood::hasRight()
123
+{
124
+    return right != NULL;
125
+}
126
+
127
+/////////////////////////////// ATOM BUCKET ////////////////////////////////////
128
+
129
+// note much of this code is unrolled intentionally so that the most frequent
130
+// cases are fast
131
+
132
+AtomBucket::AtomBucket()
133
+    : mFull(false), mOverflow(NULL), mPrev(NULL), mNext(NULL)
134
+{}
135
+
136
+unsigned AtomBucket::size() const
137
+{
138
+    unsigned thisSize = mFull ? 1 : 0;
139
+    return mOverflow == NULL ? thisSize : thisSize + mOverflow->size();
140
+}
141
+
142
+bool AtomBucket::isEmpty() const
143
+{
144
+    return !mFull;
145
+}
146
+
147
+bool AtomBucket::contains(uint64_t pos) const
148
+{
149
+    if (mOverflow != NULL && pos > mBucket.pos)
26 150
     {
27
-        pos = mRng.uniform64(0, mDomainSize);
28
-    } while (mPositionLookup.count(pos) > 0);
29
-    return pos;
151
+        return mOverflow->contains(pos);
152
+    }
153
+    else
154
+    {
155
+        return mFull ? pos == mBucket.pos : false;
156
+    }
30 157
 }
31 158
 
32
-// O(1)
33
-uint64_t AtomicDomain::size() const
159
+Atom* AtomBucket::operator[](unsigned index)
34 160
 {
35
-    return mAtoms.size();
161
+    return index == 0 ? &mBucket : mOverflow->operator[](index - 1);
36 162
 }
37 163
 
38
-// O(1)
39
-Atom& AtomicDomain::_left(const Atom &atom)
164
+void AtomBucket::insert(uint64_t pos, float mass)
40 165
 {
41
-    return mAtoms[atom.leftNdx - 1];
166
+    if (!mFull)
167
+    {
168
+        mBucket = Atom(pos, mass);
169
+        mFull = true;
170
+    }
171
+    else
172
+    {
173
+        if (mOverflow == NULL)
174
+        {
175
+            mOverflow = createAtomBucket();
176
+            mOverflow->mPrev = this;
177
+            mOverflow->mNext = mNext;
178
+        }
179
+        
180
+        if (pos < mBucket.pos)
181
+        {
182
+            mOverflow->insert(mBucket.pos, mBucket.mass);
183
+            mBucket = Atom(pos, mass);
184
+        }
185
+        else
186
+        {
187
+            mOverflow->insert(pos, mass);
188
+        }
189
+    }
42 190
 }
43 191
 
44
-// O(1)
45
-Atom& AtomicDomain::_right(const Atom &atom)
192
+// assumes pos is contained in this chain
193
+void AtomBucket::erase(uint64_t pos)
46 194
 {
47
-    return mAtoms[atom.rightNdx - 1];
195
+    GAPS_ASSERT(contains(pos));
196
+    if (mBucket.pos == pos)
197
+    {
198
+        if (mOverflow != NULL && !mOverflow->isEmpty())
199
+        {
200
+            mBucket = *(mOverflow->front());
201
+            mOverflow->eraseFront();
202
+        }
203
+        else
204
+        {
205
+            mFull = false;
206
+            connectAdjacent();
207
+        }
208
+    }
209
+    else
210
+    {
211
+        mOverflow->erase(pos);
212
+    }
48 213
 }
49 214
 
50
-// O(1)
51
-const Atom& AtomicDomain::left(const Atom &atom) const
215
+void AtomBucket::eraseFront()
52 216
 {
53
-    return mAtoms[atom.leftNdx - 1];
217
+    if (mOverflow != NULL && !mOverflow->isEmpty())
218
+    {
219
+        mBucket = *(mOverflow->front());
220
+        mOverflow->eraseFront();
221
+    }
222
+    else
223
+    {
224
+        mFull = false;
225
+        connectAdjacent();
226
+    }
54 227
 }
55 228
 
56
-// O(1)
57
-const Atom& AtomicDomain::right(const Atom &atom) const
229
+AtomNeighborhood AtomBucket::getNeighbors(unsigned index)
58 230
 {
59
-    return mAtoms[atom.rightNdx - 1];
231
+    return AtomNeighborhood(getLeft(index), this->operator[](index), getRight(index));
60 232
 }
61 233
 
62
-// O(1)
63
-bool AtomicDomain::hasLeft(const Atom &atom) const
234
+AtomNeighborhood AtomBucket::getRightNeighbor(unsigned index)
64 235
 {
65
-    return atom.leftNdx != 0;
236
+    return AtomNeighborhood(NULL, this->operator[](index), getRight(index));
66 237
 }
67 238
 
68
-// O(1)
69
-bool AtomicDomain::hasRight(const Atom &atom) const
239
+// needs to propogate through overflow
240
+void AtomBucket::setRightAdjacentBucket(AtomBucket *bucket)
70 241
 {
71
-    return atom.rightNdx != 0;
242
+    mNext = bucket;
243
+    if (mOverflow != NULL)
244
+    {
245
+        mOverflow->setRightAdjacentBucket(bucket);
246
+    }
72 247
 }
73 248
 
74
-// O(logN)
75
-Atom AtomicDomain::insert(uint64_t pos, float mass)
249
+void AtomBucket::setLeftAdjacentBucket(AtomBucket *bucket)
76 250
 {
77
-    // insert position into map
78
-    std::map<uint64_t, uint64_t>::iterator iter, iterLeft, iterRight;
79
-    iter = mAtomPositions.insert(std::pair<uint64_t, uint64_t>(pos, size())).first;
80
-    iterLeft = iter;
81
-    iterRight = iter;
251
+    mPrev = bucket;
252
+}
82 253
 
83
-    // find neighbors
84
-    Atom atom(pos, mass);
85
-    if (iter != mAtomPositions.begin())
254
+void AtomBucket::connectAdjacent()
255
+{
256
+    if (mNext != NULL)
86 257
     {
87
-        --iterLeft;
88
-        atom.leftNdx = iterLeft->second + 1;
89
-        _left(atom).rightNdx = size() + 1;
258
+        mNext->setLeftAdjacentBucket(mPrev);
90 259
     }
91
-    if (++iter != mAtomPositions.end())
260
+    if (mPrev != NULL)
92 261
     {
93
-        ++iterRight;
94
-        atom.rightNdx = iterRight->second + 1;
95
-        _right(atom).leftNdx = size() + 1;
96
-    } 
262
+        mPrev->setRightAdjacentBucket(mNext);
263
+    }
264
+}
97 265
 
98
-    // add atom to vector and lookup table
99
-    mPositionLookup.insert(std::pair<uint64_t, uint64_t>(pos, size()));
100
-    mAtoms.push_back(atom);
266
+Atom* AtomBucket::front()
267
+{
268
+    return &mBucket;
269
+}
101 270
 
102
-    return atom;
271
+Atom* AtomBucket::back()
272
+{
273
+    if (mOverflow == NULL || (mOverflow != NULL && mOverflow->isEmpty()))
274
+    {
275
+        return &mBucket;
276
+    }
277
+    return mOverflow->back();
103 278
 }
104 279
 
105
-// O(logN)
106
-// erase directly from position map and used positions hash set
107
-// swap with last atom in vector, pop off the back
108
-void AtomicDomain::erase(uint64_t pos)
280
+Atom* AtomBucket::getLeft(unsigned index)
109 281
 {
110
-    // get vector index of this atom
111
-    uint64_t index = mPositionLookup.at(pos);
282
+    if (index == 0)
283
+    {
284
+        return mPrev != NULL ? mPrev->back() : NULL;
285
+    }
286
+    else if (index == 1)
287
+    {
288
+        return &mBucket;
289
+    }
290
+    else
291
+    {
292
+        return mOverflow->getLeft(index - 1);
293
+    }    
294
+}
112 295
 
113
-    // connect neighbors of atom to be deleted
114
-    if (hasLeft(mAtoms[index]))
296
+Atom* AtomBucket::getRight(unsigned index)
297
+{
298
+    if (index == 0)
299
+    {
300
+        if (mOverflow != NULL && !mOverflow->isEmpty())
301
+        {
302
+            return mOverflow->front();
303
+        }
304
+        else
305
+        {
306
+            return mNext != NULL ? mNext->front() : NULL;
307
+        }
308
+    }
309
+    else // must be overflowed
115 310
     {
116
-        _left(mAtoms[index]).rightNdx = mAtoms[index].rightNdx;
311
+        return mOverflow->getRight(index - 1);
117 312
     }
118
-    if (hasRight(mAtoms[index]))
313
+}
314
+
315
+Archive& operator<<(Archive &ar, AtomBucket &b)
316
+{
317
+    bool hasOverflow = (b.mOverflow != NULL);
318
+    ar << b.mBucket << b.mFull << hasOverflow;
319
+
320
+    if (hasOverflow)
119 321
     {
120
-        _right(mAtoms[index]).leftNdx = mAtoms[index].leftNdx;
322
+        ar << *b.mOverflow;
121 323
     }
324
+    return ar;
325
+}
326
+
327
+Archive& operator>>(Archive& ar, AtomBucket &b)
328
+{
329
+    bool hasOverflow = false;
330
+    ar >> b.mBucket >> b.mFull >> hasOverflow;
122 331
 
123
-    // replace with atom from back
124
-    if (index < size() - 1)
332
+    if (hasOverflow)
125 333
     {
126
-        mAtoms[index] = mAtoms.back();
334
+        b.mOverflow = new AtomBucket();
335
+        ar >> *b.mOverflow;
336
+    }
337
+    return ar;
338
+}
127 339
 
128
-        // update position in map
129
-        mAtomPositions.erase(mAtoms[index].pos);
130
-        mAtomPositions.insert(std::pair<uint64_t, uint64_t>(mAtoms[index].pos,
131
-            index));
340
+////////////////////////////// ATOM HASH MAP ///////////////////////////////////
132 341
 
133
-        mPositionLookup.erase(mAtoms[index].pos);
134
-        mPositionLookup.insert(std::pair<uint64_t, uint64_t>(mAtoms[index].pos,
135
-            index));
342
+AtomHashMap::AtomHashMap(unsigned expectedNumAtoms)
343
+    :
344
+mExpectedNumAtoms(expectedNumAtoms), mLongestBucketSize(0), mSize(0),
345
+mHashMap(std::vector<AtomBucket>(expectedNumAtoms, AtomBucket()))
346
+{}
347
+
348
+void AtomHashMap::setTotalLength(uint64_t length)
349
+{
350
+    GAPS_ASSERT(length % mExpectedNumAtoms == 0);
351
+    mTotalLength = length;
352
+    mBucketLength = mTotalLength / mExpectedNumAtoms;
353
+}
354
+
355
+// pos ranges from 1 to mTotalLength
356
+unsigned AtomHashMap::hash(uint64_t pos) const
357
+{
358
+    return pos / mBucketLength;
359
+}
360
+
361
+Atom* AtomHashMap::front()
362
+{
363
+    return mHashMap[*(mFullBuckets.begin())].front();
364
+}
136 365
     
137
-        // update moved atom's neighbors
138
-        if (hasLeft(mAtoms[index]))
366
+HashMapIndex AtomHashMap::getRandomIndex() const
367
+{
368
+    GAPS_ASSERT(mSize > 0);
369
+    while (true)
370
+    {
371
+        unsigned bucket = hash(mRng.uniform64(0, mTotalLength - 1));
372
+        if (!mHashMap[bucket].isEmpty())
373
+        {
374
+            unsigned pos = mRng.uniform32(0, mLongestBucketSize - 1);
375
+            if (pos < mHashMap[bucket].size())
376
+            {
377
+                return HashMapIndex(bucket, pos);
378
+            }
379
+        }
380
+    }
381
+}
382
+
383
+Atom* AtomHashMap::randomAtom()
384
+{
385
+    HashMapIndex ndx = getRandomIndex();
386
+    return mHashMap[ndx.bucket][ndx.index];
387
+}
388
+
389
+AtomNeighborhood AtomHashMap::randomAtomWithNeighbors()
390
+{
391
+    HashMapIndex ndx = getRandomIndex();
392
+    return mHashMap[ndx.bucket].getNeighbors(ndx.index);
393
+}
394
+
395
+AtomNeighborhood AtomHashMap::randomAtomWithRightNeighbor()
396
+{
397
+    HashMapIndex ndx = getRandomIndex();
398
+    return mHashMap[ndx.bucket].getRightNeighbor(ndx.index);
399
+}
400
+
401
+bool AtomHashMap::contains(uint64_t pos) const
402
+{
403
+    unsigned index = hash(pos);
404
+    return !(mHashMap[index].isEmpty()) && mHashMap[index].contains(pos);
405
+}
406
+
407
+unsigned AtomHashMap::size() const
408
+{
409
+    return mSize;
410
+}
411
+
412
+// usually O(1), might be O(logN)
413
+void AtomHashMap::insert(uint64_t pos, float mass)
414
+{
415
+    // hash position once
416
+    unsigned index = hash(pos);
417
+
418
+    // if inserting into an empty bucket, mark bucket as non-empty
419
+    // and connect adjacent buckets O(logN)
420
+    if (mHashMap[index].isEmpty())
421
+    {
422
+        mHashMap[index].insert(pos, mass);
423
+        std::set<unsigned>::iterator it = mFullBuckets.insert(index).first;
424
+        if (it != mFullBuckets.begin())
139 425
         {
140
-            _left(mAtoms[index]).rightNdx = index + 1;
426
+            --it;
427
+            mHashMap[*it].setRightAdjacentBucket(&mHashMap[index]);
428
+            mHashMap[index].setLeftAdjacentBucket(&mHashMap[*it]);
429
+            ++it;
141 430
         }
142
-        if (hasRight(mAtoms[index]))
431
+        if (++it != mFullBuckets.end())
143 432
         {
144
-            _right(mAtoms[index]).leftNdx = index + 1;
433
+            mHashMap[*it].setLeftAdjacentBucket(&mHashMap[index]);
434
+            mHashMap[index].setRightAdjacentBucket(&mHashMap[*it]);
145 435
         }
146 436
     }
437
+    else
438
+    {
439
+        mHashMap[index].insert(pos, mass);
440
+    }
441
+    GAPS_ASSERT(mHashMap[index].contains(pos));
147 442
 
148
-    // delete atom from vector in O(1), map in O(logN)
149
-    mAtomPositions.erase(pos);
150
-    mAtoms.pop_back();
151
-    mPositionLookup.erase(pos);
443
+    // insert atom
444
+    ++mSize;
445
+
446
+    // check if this is now the longest bucket
447
+    if (mHashMap[index].size() > mLongestBucketSize)
448
+    {
449
+        mWhichLongestBucket = index;
450
+        ++mLongestBucketSize;
451
+    }
452
+}
453
+
454
+// usually O(1), sometimes O(logN), rarely O(N)
455
+void AtomHashMap::erase(uint64_t pos)
456
+{
457
+    // erase atom
458
+    unsigned index = hash(pos);
459
+    GAPS_ASSERT(!mHashMap[index].isEmpty());
460
+    GAPS_ASSERT(mHashMap[index].contains(pos));
461
+    --mSize;
462
+    
463
+    // mark bucket as empty if this was the last atom O(logN)
464
+    mHashMap[index].erase(pos);
465
+    if (mHashMap[index].isEmpty())
466
+    {
467
+        mFullBuckets.erase(index);
468
+    }
469
+    GAPS_ASSERT(!mHashMap[index].contains(pos));
470
+
471
+    // if this atom was in the largest bucket, find the new largest bucket O(N)
472
+    if (index == mWhichLongestBucket)
473
+    {
474
+        --mLongestBucketSize;
475
+        std::set<unsigned>::iterator it = mFullBuckets.begin();
476
+        for (; it != mFullBuckets.end(); ++it)
477
+        {
478
+            if (mHashMap[*it].size() > mLongestBucketSize)
479
+            {
480
+                mLongestBucketSize = mHashMap[*it].size();
481
+                mWhichLongestBucket = *it;
482
+            }
483
+        }
484
+    }
485
+}
486
+
487
+void AtomHashMap::move(uint64_t src, uint64_t dest, float mass)
488
+{
489
+    unsigned srcHash = hash(src);
490
+    unsigned destHash = hash(dest);
491
+    if (srcHash != destHash)
492
+    {
493
+        erase(src);
494
+        insert(dest, mass);
495
+    }
496
+}
497
+
498
+Archive& operator<<(Archive& ar, AtomHashMap &h)
499
+{
500
+    ar << h.mExpectedNumAtoms << h.mBucketLength << h.mTotalLength <<
501
+        h.mWhichLongestBucket << h.mLongestBucketSize << h.mSize;
502
+
503
+    ar << h.mFullBuckets.size();
504
+    std::set<unsigned>::iterator it = h.mFullBuckets.begin();
505
+    for (; it != h.mFullBuckets.end(); ++it)
506
+    {
507
+        ar << *it << h.mHashMap[*it];
508
+    }
509
+    return ar;
510
+}
511
+
512
+Archive& operator>>(Archive& ar, AtomHashMap &h)
513
+{
514
+    ar >> h.mExpectedNumAtoms >> h.mBucketLength >> h.mTotalLength >>
515
+        h.mWhichLongestBucket >> h.mLongestBucketSize >> h.mSize;
516
+
517
+    unsigned nBuckets = 0;
518
+    ar >> nBuckets;
519
+
520
+    for (unsigned i = 0; i < nBuckets; ++i)
521
+    {
522
+        unsigned thisBucket = 0;
523
+        ar >> thisBucket;
524
+        h.mFullBuckets.insert(thisBucket);
525
+        ar >> h.mHashMap[thisBucket];
526
+    }
527
+    return ar;
528
+}
529
+
530
+////////////////////////////// ATOMIC DOMAIN ///////////////////////////////////
531
+
532
+AtomicDomain::AtomicDomain(unsigned nBins) : mAtoms(nBins)
533
+{
534
+    uint64_t binLength = std::numeric_limits<uint64_t>::max() / nBins;
535
+    mDomainLength = binLength * nBins;
536
+    mAtoms.setTotalLength(mDomainLength);
537
+}
538
+
539
+void AtomicDomain::setDomainSize(uint64_t size)
540
+{
541
+    // do nothing
542
+}
543
+
544
+Atom* AtomicDomain::front()
545
+{
546
+    return mAtoms.front();
547
+}
548
+
549
+Atom* AtomicDomain::randomAtom()
550
+{
551
+    return mAtoms.randomAtom();
552
+}
553
+
554
+AtomNeighborhood AtomicDomain::randomAtomWithNeighbors()
555
+{
556
+    return mAtoms.randomAtomWithNeighbors();
557
+}
558
+
559
+AtomNeighborhood AtomicDomain::randomAtomWithRightNeighbor()
560
+{
561
+    return mAtoms.randomAtomWithRightNeighbor();
562
+}
563
+
564
+uint64_t AtomicDomain::randomFreePosition() const
565
+{
566
+    uint64_t pos = mRng.uniform64(1, mDomainLength);
567
+    while (mAtoms.contains(pos))
568
+    {
569
+        pos = mRng.uniform64(1, mDomainLength);
570
+    } 
571
+    return pos;
572
+}
573
+
574
+uint64_t AtomicDomain::size() const
575
+{
576
+    return mAtoms.size();
152 577
 }
153 578
 
154 579
 void AtomicDomain::cacheInsert(uint64_t pos, float mass) const
... ...
@@ -158,7 +583,7 @@ void AtomicDomain::cacheInsert(uint64_t pos, float mass) const
158 583
     {
159 584
         ndx = mInsertCacheIndex++;
160 585
     }
161
-    mInsertCache[ndx] = RawAtom(pos, mass);
586
+    mInsertCache[ndx] = Atom(pos, mass);
162 587
 }
163 588
 
164 589
 void AtomicDomain::cacheErase(uint64_t pos) const
... ...
@@ -171,70 +596,54 @@ void AtomicDomain::cacheErase(uint64_t pos) const
171 596
     mEraseCache[ndx] = pos;
172 597
 }
173 598
 
599
+void AtomicDomain::cacheMove(uint64_t src, uint64_t dest, float mass) const
600
+{
601
+    unsigned ndx = 0;
602
+    #pragma omp critical(atomicMove)
603
+    {
604
+        ndx = mMoveCacheIndex++;
605
+    }
606
+    mMoveCache[ndx] = MoveOperation(src, dest, mass);
607
+}
608
+
174 609
 void AtomicDomain::resetCache(unsigned n)
175 610
 {
176 611
     mInsertCacheIndex = 0;
177 612
     mEraseCacheIndex = 0;
613
+    mMoveCacheIndex = 0;
178 614
     mInsertCache.resize(n);
179 615
     mEraseCache.resize(n);
616
+    mMoveCache.resize(n);
180 617
 }
181 618
 
182 619
 void AtomicDomain::flushCache()
183 620
 {
184 621
     for (unsigned i = 0; i < mEraseCacheIndex; ++i)
185 622
     {
186
-        erase(mEraseCache[i]);
623
+        mAtoms.erase(mEraseCache[i]);
187 624
     }
188 625
 
189 626
     for (unsigned i = 0; i < mInsertCacheIndex; ++i)
190 627
     {
191
-        insert(mInsertCache[i].pos, mInsertCache[i].mass);
628
+        mAtoms.insert(mInsertCache[i].pos, mInsertCache[i].mass);
629
+    }
630
+
631
+    for (unsigned i = 0; i < mMoveCacheIndex; ++i)
632
+    {
633
+        mAtoms.move(mMoveCache[i].src, mMoveCache[i].dest, mMoveCache[i].mass);
192 634
     }
193 635
 
194 636
     mInsertCache.clear();
195 637
     mEraseCache.clear();
196
-}
197
-
198
-// O(1)
199
-void AtomicDomain::updateMass(uint64_t pos, float newMass)
200
-{
201
-    mAtoms[mPositionLookup.at(pos)].mass = newMass;
202
-}
203
-
204
-Archive& operator<<(Archive &ar, Atom &a)
205
-{
206
-    ar << a.leftNdx << a.rightNdx << a.pos << a.mass;
207
-    return ar;
208
-}
209
-
210
-Archive& operator>>(Archive &ar, Atom &a)
211
-{
212
-    ar >> a.leftNdx >> a.rightNdx >> a.pos >> a.mass;
213
-    return ar;
638
+    mMoveCache.clear();
214 639
 }
215 640
 
216 641
 Archive& operator<<(Archive &ar, AtomicDomain &domain)
217 642
 {
218
-    unsigned nAtoms = domain.mAtoms.size();
219
-    ar << nAtoms << domain.mDomainSize << domain.mRng;
220
-
221
-    for (unsigned i = 0; i < nAtoms; ++i)
222
-    {
223
-        ar << domain.mAtoms[i];
224
-    }
225
-    return ar;
643
+    ar << domain.mDomainLength << domain.mAtoms << domain.mRng;
226 644
 }
227 645
 
228 646
 Archive& operator>>(Archive &ar, AtomicDomain &domain)
229 647
 {
230
-    unsigned nAtoms = 0;
231
-    ar >> nAtoms >> domain.mDomainSize >> domain.mRng;
232
-
233
-    Atom a(0, 0.f);
234
-    for (unsigned i = 0; i < nAtoms; ++i)
235
-    {
236
-        ar >> a;
237
-        domain.insert(a.pos, a.mass);
238
-    }
239
-    return ar;
648
+    ar >> domain.mDomainLength >> domain.mAtoms << domain.mRng;
240 649
 }
... ...
@@ -4,109 +4,218 @@
4 4
 #include "Archive.h"
5 5
 #include "math/Random.h"
6 6
 
7
-#include <boost/unordered_map.hpp>
8
-
9
-#include <cstddef>
10
-#include <map>
11
-#include <stdint.h>
12 7
 #include <vector>
8
+#include <set>
13 9
 
14
-class AtomicDomain;
15
-
16
-// Want the neighbors to be pointers, but can't use pointers since vector
17
-// is resized, shifted integers allow for 0 to be "null" and 1 to be the
18
-// first index. AtomicDomain is responsible for managing this correctly.
19
-// TODO use get/set neighbor
20 10
 struct Atom
21 11
 {
22
-private:    
12
+    uint64_t pos;
13
+    float mass;
23 14
 
24
-    friend class AtomicDomain;
25
-    uint64_t leftNdx; // shift these up 1, 0 == no neighbor
26
-    uint64_t rightNdx;
15
+    Atom();
16
+    Atom(uint64_t p, float m);
27 17
 
28
-public:    
18
+    void operator=(Atom other);
29 19
 
30
-    uint64_t pos;
31
-    float mass;
20
+    friend Archive& operator<<(Archive& ar, Atom &a);
21
+    friend Archive& operator>>(Archive& ar, Atom &a);
22
+};
32 23
 
33
-    Atom(uint64_t p, float m)
34
-        : leftNdx(0), rightNdx(0), pos(p), mass(m)
35
-    {}
24
+struct AtomNeighborhood
25
+{
26
+    Atom* center;
27
+    Atom* left;
28
+    Atom* right;
36 29
 
37
-    bool operator==(const Atom &other) const
38
-    {
39
-        return pos == other.pos;
40
-    }
30
+    AtomNeighborhood(Atom *l, Atom *c, Atom *r);
41 31
 
42
-    bool operator!=(const Atom &other) const
43
-    {
44
-        return !(pos == other.pos);
45
-    }
32
+    bool hasLeft();
33
+    bool hasRight();
34
+};
46 35
 
47
-    // serialization
48
-    friend Archive& operator<<(Archive &ar, Atom &a);
49
-    friend Archive& operator>>(Archive &ar, Atom &a);
36
+class AtomBucket
37
+{
38
+public:
39
+
40
+    AtomBucket();
41
+
42
+    unsigned size() const;
43
+    bool isEmpty() const;
44
+    bool contains(uint64_t pos) const;
45
+
46
+    Atom* front();
47
+    Atom* operator[](unsigned index);
48
+
49
+    void insert(uint64_t pos, float mass);
50
+    void erase(uint64_t pos);
51
+        
52
+    AtomNeighborhood getNeighbors(unsigned index);
53
+    AtomNeighborhood getRightNeighbor(unsigned index);
54
+
55
+    void setRightAdjacentBucket(AtomBucket *bucket);
56
+    void setLeftAdjacentBucket(AtomBucket *bucket);
57
+
58
+    AtomBucket& operator=(const AtomBucket& other);
59
+
60
+    unsigned getIndex(uint64_t pos);
61
+    Atom* getLeft(unsigned index);
62
+    Atom* getRight(unsigned index);
63
+
64
+#ifndef GAPS_INTERNAL_TESTS
65
+private:
66
+#endif
67
+
68
+    Atom mBucket;
69
+    bool mFull;
70
+    
71
+    AtomBucket *mOverflow;
72
+    AtomBucket *mPrev;
73
+    AtomBucket *mNext;
74
+
75
+    void eraseFront();
76
+    void connectAdjacent();
77
+
78
+    Atom* back();
79
+    
80
+    friend Archive& operator<<(Archive& ar, AtomBucket &b);
81
+    friend Archive& operator>>(Archive& ar, AtomBucket &b);
50 82
 };
51 83
 
52
-struct RawAtom
84
+struct HashMapIndex
53 85
 {
54
-    uint64_t pos;
55
-    float mass;
86
+    unsigned bucket;
87
+    unsigned index;
56 88
 
57
-    RawAtom() : pos(0), mass(0.f) {}
58
-    RawAtom(uint64_t p, float m) : pos(p), mass(m) {}
89
+    HashMapIndex(unsigned b, unsigned i) : bucket(b), index(i) {}
59 90
 };
60 91
 
92
+// hash map of chained AtomBuckets
61 93
 // data structure that holds atoms
62
-class AtomicDomain
94
+// accessing happens much more than insert/erase so more time is spent
95
+// up front (sorting) to save time on access
96
+// note that atoms can never occupy position 0
97
+class AtomHashMap
63 98
 {
99
+public:
100
+
101
+    // required that length is a multiple of expectedNumAtoms
102
+    AtomHashMap(unsigned expectedNumAtoms);
103
+    void setTotalLength(uint64_t length);
104
+
105
+    // TODO while moving the atoms cannot change the atom order, it can
106
+    // move an atom from one bucket to another - this moved atom must be
107
+    // the front atom moving to the back of the prev bucket or the back
108
+    // atom moving to the front of the next atom - effectively creating
109
+    // a reorder event - it may even need allocation
110
+
111
+    Atom* front();
112
+    Atom* randomAtom();
113
+    AtomNeighborhood randomAtomWithNeighbors();
114
+    AtomNeighborhood randomAtomWithRightNeighbor();
115
+
116
+    bool contains(uint64_t pos) const;
117
+    unsigned size() const;
118
+
119
+    void insert(uint64_t pos, float mass);
120
+    void erase(uint64_t pos);
121
+    void move(uint64_t src, uint64_t dest, float mass); // cannot reorder
122
+
123
+    Atom* getLeft(uint64_t pos);
124
+    Atom* getRight(uint64_t pos);
125
+
126
+    bool hasLeft(uint64_t pos);
127
+    bool hasRight(uint64_t pos);
128
+
129
+    void updateMass(uint64_t pos, float newMass);
130
+
131
+#ifndef GAPS_INTERNAL_TESTS
64 132
 private:
133
+#endif
65 134
 
66
-    // size of atomic domain
67
-    uint64_t mDomainSize;
135
+    unsigned mSize;
136
+    unsigned mWhichLongestBucket;
137
+    unsigned mLongestBucketSize;
68 138
 
69
-    // domain storage
70
-    std::vector<Atom> mAtoms;
71
-    std::map<uint64_t, uint64_t> mAtomPositions;
72
-    boost::unordered_map<uint64_t, uint64_t> mPositionLookup;
139
+    unsigned mExpectedNumAtoms;
73 140
 
74
-    mutable std::vector<RawAtom> mInsertCache;
75
-    mutable std::vector<uint64_t> mEraseCache;
141
+    uint64_t mBucketLength;
142
+    uint64_t mTotalLength;
76 143
 
77
-    mutable unsigned mInsertCacheIndex;
78
-    mutable unsigned mEraseCacheIndex;
144
+    std::vector<AtomBucket> mHashMap;
145
+    std::set<unsigned> mFullBuckets; // TODO use IntFixedHashSet
79 146
 
147
+    // random number generator
80 148
     mutable GapsRng mRng;
81 149
 
82
-    Atom& _left(const Atom &atom);
83
-    Atom& _right(const Atom &atom);
150
+    unsigned hash(uint64_t pos) const;
151
+    HashMapIndex getRandomIndex() const;
84 152
 
85
-    void erase(uint64_t pos);
86
-    Atom insert(uint64_t pos, float mass);
153
+    friend Archive& operator<<(Archive& ar, AtomHashMap &h);
154
+    friend Archive& operator>>(Archive& ar, AtomHashMap &h);
155
+};
156
+
157
+struct MoveOperation
158
+{
159
+    uint64_t src;
160
+    uint64_t dest;
161
+    float mass;
162
+
163
+    MoveOperation() : src(0), dest(0), mass(0.f) {}
164
+
165
+    MoveOperation(uint64_t s, uint64_t d, float m) :
166
+        src(s), dest(d), mass(m)
167
+    {}
168
+};
87 169
 
170
+class AtomicDomain
171
+{
88 172
 public:
89 173
 
90
-    void setDomainSize(uint64_t size) { mDomainSize = size; }
174
+    AtomicDomain(unsigned nBins);
175
+
176
+    void setDomainSize(uint64_t size);
91 177
 
92 178
     // access atoms
93
-    Atom front() const;
94
-    Atom randomAtom() const;
179
+    Atom* front();
180
+    Atom* randomAtom();
181
+    AtomNeighborhood randomAtomWithNeighbors();
182
+    AtomNeighborhood randomAtomWithRightNeighbor();
183
+
95 184
     uint64_t randomFreePosition() const;
96 185
     uint64_t size() const;
97 186
 
98
-    const Atom& left(const Atom &atom) const;
99
-    const Atom& right(const Atom &atom) const;
100
-    bool hasLeft(const Atom &atom) const;
101
-    bool hasRight(const Atom &atom) const;
102
-
103 187
     // modify domain
104 188
     void cacheInsert(uint64_t pos, float mass) const;
105 189
     void cacheErase(uint64_t pos) const;
106
-    void updateMass(uint64_t pos, float newMass);
190
+    void cacheMove(uint64_t src, uint64_t dest, float mass) const;
107 191
     void flushCache();
108 192
     void resetCache(unsigned n);
109 193
 
194
+    // serialization
195
+    friend Archive& operator<<(Archive &ar, AtomicDomain &domain);
196
+    friend Archive& operator>>(Archive &ar, AtomicDomain &domain);
197
+
198
+private:
199
+
200
+    // size of atomic domain to ensure all bins are equal length
201
+    uint64_t mDomainLength;
202
+
203
+    // domain storage, specialized hash map
204
+    mutable AtomHashMap mAtoms;
205
+
206
+    // holds cache of operations
207
+    mutable std::vector<Atom> mInsertCache;
208
+    mutable std::vector<uint64_t> mEraseCache;
209
+    mutable std::vector<MoveOperation> mMoveCache;
210
+
211
+    // current index in the operation cache
212
+    mutable unsigned mInsertCacheIndex;
213
+    mutable unsigned mEraseCacheIndex;
214
+    mutable unsigned mMoveCacheIndex;
215
+
216
+    // random number generator
217
+    mutable GapsRng mRng;
218
+
110 219
     // serialization
111 220
     friend Archive& operator<<(Archive &ar, AtomicDomain &domain);
112 221
     friend Archive& operator>>(Archive &ar, AtomicDomain &domain);
... ...
@@ -60,21 +60,19 @@ protected:
60 60
 
61 61
     T* impl();
62 62
 
63
-    void processProposal(const AtomicProposal &prop);
63
+    void processProposal(AtomicProposal *prop);
64 64
 
65 65
     void addMass(uint64_t pos, float mass, unsigned row, unsigned col);
66 66
     void removeMass(uint64_t pos, float mass, unsigned row, unsigned col);
67 67
 
68
-    void birth(uint64_t pos, unsigned row, unsigned col, GapsRng *rng);
69
-    void death(uint64_t pos, float mass, unsigned row, unsigned col, GapsRng *rng);
70
-    void move(uint64_t src, float mass, uint64_t dest, unsigned r1, unsigned c1,
71
-        unsigned r2, unsigned c2, GapsRng *rng);
72
-    void exchange(uint64_t p1, float m1, uint64_t p2, float m2, unsigned r1,
73
-        unsigned c1, unsigned r2, unsigned c2, GapsRng *rng);
68
+    void birth(AtomicProposal *prop);
69
+    void death(AtomicProposal *prop);
70
+    void move(AtomicProposal *prop);
71
+    void exchange(AtomicProposal *prop);
74 72
 
75
-    bool updateAtomMass(uint64_t pos, float mass, float delta);
76
-    void acceptExchange(uint64_t p1, float m1, float d1, uint64_t p2, float m2,
77
-        float d2, unsigned r1, unsigned c1, unsigned r2, unsigned c2);
73
+    bool updateAtomMass(Atom *atom, float delta);
74
+    void acceptExchange(AtomicProposal *prop, float d1, unsigned r1,
75
+        unsigned c1, unsigned r2, unsigned c2);
78 76
 
79 77
     std::pair<float, bool> gibbsMass(AlphaParameters alpha, GapsRng *rng);
80 78
     std::pair<float, bool> gibbsMass(AlphaParameters alpha, float m1, float m2,
... ...
@@ -210,7 +208,8 @@ mDMatrix(data, transposeData, partitionRows, indices),
210 208
 mSMatrix(mDMatrix.pmax(0.1f, 0.1f)),
211 209
 mAPMatrix(mDMatrix.nRow(), mDMatrix.nCol()),
212 210
 mMatrix(amp ? mDMatrix.nRow() : nPatterns, amp ? nPatterns : mDMatrix.nCol()),
213
-mOtherMatrix(NULL), mQueue(mMatrix.nRow() * mMatrix.nCol()), mLambda(0.f),
211
+mOtherMatrix(NULL), mQueue(mMatrix.nRow() * mMatrix.nCol()),
212
+mDomain(mMatrix.nRow() * mMatrix.nCol()), mLambda(0.f),
214 213
 mMaxGibbsMass(100.f), mAnnealingTemp(1.f), mNumRows(mMatrix.nRow()),
215 214
 mNumCols(mMatrix.nCol()), mAvgQueue(0.f), mNumQueues(0.f)
216 215
 {
... ...
@@ -286,7 +285,7 @@ void GibbsSampler<T, MatA, MatB>::update(unsigned nSteps, unsigned nCores)
286 285
         #pragma omp parallel for num_threads(nCores)
287 286
         for (unsigned i = 0; i < mQueue.size(); ++i)
288 287
         {
289
-            processProposal(mQueue[i]);
288
+            processProposal(&mQueue[i]);
290 289
         }
291 290
         mDomain.flushCache();
292 291
         mQueue.clear();
... ...
@@ -327,41 +326,7 @@ float GibbsSampler<T, MatA, MatB>::getAvgQueue() const
327 326
 template <class T, class MatA, class MatB>
328 327
 bool GibbsSampler<T, MatA, MatB>::internallyConsistent()
329 328
 {
330
-    if (mDomain.size() > 0)
331
-    {
332
-        Atom a = mDomain.front();
333
-        float current = a.mass;
334
-        uint64_t row = impl()->getRow(a.pos);
335
-        uint64_t col = impl()->getCol(a.pos);
336
-
337
-        while (mDomain.hasRight(a))
338
-        {
339
-            a = mDomain.right(a);
340
-            if (row != impl()->getRow(a.pos) || col != impl()->getCol(a.pos))
341
-            {
342
-                float matVal = mMatrix(row, col);
343
-                if (std::abs(current - matVal) > 0.1f)
344
-                {
345
-                    gaps_printf("mass difference detected at row %lu, column %lu: %f %f\n",
346
-                        row, col, current, matVal); 
347
-                    return false;
348
-                }
349
-                
350
-                row = impl()->getRow(a.pos);
351
-                col = impl()->getCol(a.pos);
352
-                current = a.mass;
353
-            }
354
-            else
355
-            {
356
-                current += a.mass;
357
-            }
358
-        }
359
-        return true;
360
-    }
361
-    else
362
-    {
363
-        return gaps::algo::sum(mMatrix) == 0.f;
364
-    }
329
+    return true;
365 330
 }
366 331
 #endif
367 332
 
... ...
@@ -392,24 +357,22 @@ T* GibbsSampler<T, MatA, MatB>::impl()
392 357
 }
393 358
 
394 359
 template <class T, class MatA, class MatB>
395
-void GibbsSampler<T, MatA, MatB>::processProposal(const AtomicProposal &prop)
360
+void GibbsSampler<T, MatA, MatB>::processProposal(AtomicProposal *prop)
396 361
 {
397
-    unsigned r1 = impl()->getRow(prop.pos1), c1 = impl()->getCol(prop.pos1);
398
-    unsigned r2 = impl()->getRow(prop.pos2), c2 = impl()->getCol(prop.pos2);
399
-
400
-    switch (prop.type)
362
+    unsigned r1 = 0, c1 = 0, r2 = 0, c2 = 0;
363
+    switch (prop->type)
401 364
     {
402 365
         case 'B':
403
-            birth(prop.pos1, r1, c1, &prop.rng);
366
+            birth(prop);
404 367
             break;
405 368
         case 'D':
406
-            death(prop.pos1, prop.mass1, r1, c1, &prop.rng);
369
+            death(prop);
407 370
             break;
408 371
         case 'M':
409
-            move(prop.pos1, prop.mass1, prop.pos2, r1, c1, r2, c2, &prop.rng);
372
+            move(prop);
410 373
             break;
411 374
         case 'E':
412
-            exchange(prop.pos1, prop.mass1, prop.pos2, prop.mass2, r1, c1, r2, c2, &prop.rng);
375
+            exchange(prop);
413 376
             break;
414 377
     }
415 378
 }
... ...
@@ -433,25 +396,29 @@ void GibbsSampler<T, MatA, MatB>::removeMass(uint64_t pos, float mass, unsigned
433 396
 // add an atom at pos, calculate mass either with an exponential distribution
434 397
 // or with the gibbs mass calculation
435 398
 template <class T, class MatA, class MatB>
436
-void GibbsSampler<T, MatA, MatB>::birth(uint64_t pos, unsigned row,
437
-unsigned col, GapsRng *rng)
399
+void GibbsSampler<T, MatA, MatB>::birth(AtomicProposal *prop)
438 400
 {
401
+    unsigned row = impl()->getRow(prop->birthPos);
402
+    unsigned col = impl()->getCol(prop->birthPos);
403
+
439 404
     // calculate proposed mass
440 405
     float mass = 0.f;
441 406
     if (impl()->canUseGibbs(row, col))
442 407
     {
443 408
         AlphaParameters alpha = impl()->alphaParameters(row, col);
444
-        mass = gibbsMass(alpha, rng).first;
409
+        mass = gibbsMass(alpha, &(prop->rng)).first;
445 410
     }
446 411
     else
447 412
     {
448
-        mass = rng->exponential(mLambda);
413
+        mass = prop->rng.exponential(mLambda);
449 414
     }
450 415
 
451 416
     // accept mass as long as it's non-zero
452 417
     if (mass >= gaps::epsilon)
453 418
     {
454
-        addMass(pos, mass, row, col);
419
+        mDomain.cacheInsert(prop->birthPos, mass);
420
+        mMatrix(row, col) += mass;
421
+        impl()->updateAPMatrix(row, col, mass);
455 422
         mQueue.acceptBirth();
456 423
     }
457 424
     else
... ...
@@ -463,19 +430,22 @@ unsigned col, GapsRng *rng)
463 430
 // automatically accept death, attempt a rebirth at the same position, using
464 431
 // the original mass or the gibbs mass calculation
465 432
 template <class T, class MatA, class MatB>
466
-void GibbsSampler<T, MatA, MatB>::death(uint64_t pos, float mass, unsigned row,
467
-unsigned col, GapsRng *rng)
433
+void GibbsSampler<T, MatA, MatB>::death(AtomicProposal *prop)
468 434
 {
435
+    // calculate bin for this atom
436
+    unsigned row = impl()->getRow(prop->atom1->pos);
437
+    unsigned col = impl()->getCol(prop->atom1->pos);
438
+
469 439
     // kill off atom
470
-    mMatrix(row, col) = gaps::max(mMatrix(row, col) - mass, 0.f);
471
-    impl()->updateAPMatrix(row, col, -mass);
440
+    mMatrix(row, col) = gaps::max(mMatrix(row, col) - prop->atom1->mass, 0.f);
441
+    impl()->updateAPMatrix(row, col, -1.f * prop->atom1->mass);
472 442
 
473 443
     // calculate rebirth mass
474
-    float rebirthMass = mass;
444
+    float rebirthMass = prop->atom1->mass;
475 445
     AlphaParameters alpha = impl()->alphaParameters(row, col);
476 446
     if (impl()->canUseGibbs(row, col))
477 447
     {
478
-        std::pair<float, bool> gMass = gibbsMass(alpha, rng);
448
+        std::pair<float, bool> gMass = gibbsMass(alpha, &(prop->rng));
479 449
         if (gMass.second)
480 450
         {
481 451
             rebirthMass = gMass.first;
... ...
@@ -484,32 +454,42 @@ unsigned col, GapsRng *rng)
484 454
 
485 455
     // accept/reject rebirth
486 456
     float deltaLL = rebirthMass * (alpha.su - alpha.s * rebirthMass / 2.f);
487
-    if (deltaLL * mAnnealingTemp >= std::log(rng->uniform()))
457
+    if (deltaLL * mAnnealingTemp >= std::log(prop->rng.uniform()))
488 458
     {
489
-        mDomain.updateMass(pos, rebirthMass);
459
+        prop->atom1->mass = rebirthMass;
490 460
         mMatrix(row, col) += rebirthMass;
491 461
         impl()->updateAPMatrix(row, col, rebirthMass);
492 462
         mQueue.rejectDeath();
493 463
     }
494 464
     else
495 465
     {
496
-        mDomain.cacheErase(pos);
466
+        mDomain.cacheErase(prop->atom1->pos);
497 467
         mQueue.acceptDeath();
498 468
     }
499 469
 }
500 470
 
501 471
 // move mass from src to dest in the atomic domain
502 472
 template <class T, class MatA, class MatB>
503
-void GibbsSampler<T, MatA, MatB>::move(uint64_t src, float mass, uint64_t dest,
504
-unsigned r1, unsigned c1, unsigned r2, unsigned c2, GapsRng *rng)
473
+void GibbsSampler<T, MatA, MatB>::move(AtomicProposal *prop)
505 474
 {
475
+    unsigned r1 = impl()->getRow(prop->atom1->pos);
476
+    unsigned c1 = impl()->getCol(prop->atom1->pos);
477
+    unsigned r2 = impl()->getRow(prop->moveDest);
478
+    unsigned c2 = impl()->getCol(prop->moveDest);
479
+
506 480
     if (r1 != r2 || c1 != c2) // automatically reject if change in same bin
507 481
     {
508
-        float deltaLL = impl()->computeDeltaLL(r1, c1, -mass, r2, c2, mass);
509
-        if (deltaLL * mAnnealingTemp > std::log(rng->uniform()))
482
+        float deltaLL = impl()->computeDeltaLL(r1, c1, -1.f * prop->atom1->mass,
483
+            r2, c2, prop->atom1->mass);
484
+        if (deltaLL * mAnnealingTemp > std::log(prop->rng.uniform()))
510 485
         {
511
-            removeMass(src, mass, r1, c1);
512
-            addMass(dest, mass, r2, c2);
486
+            mDomain.cacheMove(prop->atom1->pos, prop->moveDest, prop->atom1->mass);
487
+
488
+            mMatrix(r1, c1) = gaps::max(mMatrix(r1, c1) - prop->atom1->mass, 0.f);
489
+            mMatrix(r2, c2) += prop->atom1->mass;
490
+
491
+            impl()->updateAPMatrix(r1, c1, -1.f * prop->atom1->mass);
492
+            impl()->updateAPMatrix(r2, c2, prop->atom1->mass);
513 493
         }
514 494
     }
515 495
 }
... ...
@@ -518,25 +498,31 @@ unsigned r1, unsigned c1, unsigned r2, unsigned c2, GapsRng *rng)
518 498
 // for one of the atoms to be deleted if it's mass becomes too small after
519 499
 // the exchange
520 500
 template <class T, class MatA, class MatB>
521
-void GibbsSampler<T, MatA, MatB>::exchange(uint64_t p1, float m1, uint64_t p2,
522
-float m2, unsigned r1, unsigned c1, unsigned r2, unsigned c2, GapsRng *rng)
501
+void GibbsSampler<T, MatA, MatB>::exchange(AtomicProposal *prop)
523 502
 {
503
+    unsigned r1 = impl()->getRow(prop->atom1->pos);
504
+    unsigned c1 = impl()->getCol(prop->atom1->pos);
505
+    unsigned r2 = impl()->getRow(prop->atom2->pos);
506
+    unsigned c2 = impl()->getCol(prop->atom2->pos);
507
+
524 508
     if (r1 != r2 || c1 != c2) // automatically reject if change in same bin
525 509
     {
510
+        float m1 = prop->atom1->mass;
511
+        float m2 = prop->atom2->mass;
512
+
526 513
         if (impl()->canUseGibbs(r1, c1, r2, c2))
527 514
         {
528 515
             AlphaParameters alpha = impl()->alphaParameters(r1, c1, r2, c2);
529
-            std::pair<float, bool> gMass = gibbsMass(alpha, m1, m2, rng);
516
+            std::pair<float, bool> gMass = gibbsMass(alpha, m1, m2, &(prop->rng));
530 517
             if (gMass.second)
531 518
             {
532
-                acceptExchange(p1, m1, gMass.first, p2, m2, -gMass.first, r1,
533
-                    c1, r2, c2);
519
+                acceptExchange(prop, gMass.first, r1, c1, r2, c2);
534 520
                 return;
535 521
             }
536 522
         }
537 523
 
538 524
         float pUpper = gaps::p_gamma(m1 + m2, 2.f, 1.f / mLambda);
539
-        float newMass = rng->inverseGammaSample(0.f, pUpper, 2.f, 1.f / mLambda);
525
+        float newMass = prop->rng.inverseGammaSample(0.f, pUpper, 2.f, 1.f / mLambda);
540 526
 
541 527
         float delta = m1 > m2 ? newMass - m1 : m2 - newMass; // change larger mass
542 528
         float pOldMass = 2.f * newMass > m1 + m2 ? gaps::max(m1, m2) : gaps::min(m1, m2);
... ...
@@ -546,15 +532,15 @@ float m2, unsigned r1, unsigned c1, unsigned r2, unsigned c2, GapsRng *rng)
546 532
 
547 533
         if (pOld == 0.f && pNew != 0.f) // special case
548 534
         {
549
-            acceptExchange(p1, m1, delta, p2, m2, -delta, r1, c1, r2, c2);
535
+            acceptExchange(prop, delta, r1, c1, r2, c2);
550 536
             return;
551 537
         }
552 538
         float deltaLL = impl()->computeDeltaLL(r1, c1, delta, r2, c2, -delta);
553 539
         float priorLL = (pOld == 0.f) ? 1.f : pOld / pNew;
554
-        float u = std::log(rng->uniform() * priorLL);
540
+        float u = std::log(prop->rng.uniform() * priorLL);
555 541
         if (u < deltaLL * mAnnealingTemp)
556 542
         {
557
-            acceptExchange(p1, m1, delta, p2, m2, -delta, r1, c1, r2, c2);
543
+            acceptExchange(prop, delta, r1, c1, r2, c2);
558 544
             return;
559 545
         }
560 546
     }
... ...
@@ -564,33 +550,32 @@ float m2, unsigned r1, unsigned c1, unsigned r2, unsigned c2, GapsRng *rng)
564 550
 // helper function for acceptExchange, used to conditionally update the mass
565 551
 // at a single position, deleting the atom if the new mass is too small
566 552
 template <class T, class MatA, class MatB>
567
-bool GibbsSampler<T, MatA, MatB>::updateAtomMass(uint64_t pos, float mass,
568
-float delta)
553
+bool GibbsSampler<T, MatA, MatB>::updateAtomMass(Atom *atom, float delta)
569 554
 {
570
-    if (mass + delta < gaps::epsilon)
555
+    if (atom->mass + delta < gaps::epsilon)
571 556
     {
572
-        mDomain.cacheErase(pos);
557
+        mDomain.cacheErase(atom->pos);
573 558
         mQueue.acceptDeath();
574 559
         return false;
575 560
     }
576
-    mDomain.updateMass(pos, mass + delta);
561
+    atom->mass += delta;
577 562
     return true;
578 563
 }
579 564
 
580 565
 // helper function for exchange step, updates the atomic domain, matrix, and
581 566
 // A*P matrix
582 567
 template <class T, class MatA, class MatB>
583
-void GibbsSampler<T, MatA, MatB>::acceptExchange(uint64_t p1, float m1,
584
-float d1, uint64_t p2, float m2, float d2, unsigned r1, unsigned c1,
585
-unsigned r2, unsigned c2)
568
+void GibbsSampler<T, MatA, MatB>::acceptExchange(AtomicProposal *prop,
569
+float d1, unsigned r1, unsigned c1, unsigned r2, unsigned c2)
586 570
 {
587
-    bool b1 = updateAtomMass(p1, m1, d1);
588
-    bool b2 = updateAtomMass(p2, m2, d2);
571
+    float d2 = -1.f * d1;
572
+    bool b1 = updateAtomMass(prop->atom1, d1);
573
+    bool b2 = updateAtomMass(prop->atom2, d2);
589 574
     GAPS_ASSERT(b1 || b2);
590 575
     
591 576
     // delete entire atom if resize would make it too small
592
-    if (!b1) { d1 = -m1; }
593
-    if (!b2) { d2 = -m2; }
577
+    if (!b1) { d1 = -1.f * prop->atom1->mass; }
578
+    if (!b2) { d2 = -1.f * prop->atom2->mass; }
594 579
 
595 580
     if (b1 && b2)
596 581
     {
... ...
@@ -598,9 +583,7 @@ unsigned r2, unsigned c2)
598 583
     }
599 584
 
600 585
     mMatrix(r1, c1) += d1;
601
-    GAPS_ASSERT(mMatrix(r1, c1) >= 0);
602 586
     mMatrix(r2, c2) += d2;
603
-    GAPS_ASSERT(mMatrix(r2, c2) >= 0);
604 587
     impl()->updateAPMatrix(r1, c1, d1);
605 588
     impl()->updateAPMatrix(r2, c2, d2);
606 589
 }
... ...
@@ -50,7 +50,7 @@ unsigned ProposalQueue::size() const
50 50
     return mQueue.size();
51 51
 }
52 52
 
53
-const AtomicProposal& ProposalQueue::operator[](int n) const
53
+AtomicProposal& ProposalQueue::operator[](int n)
54 54
 {
55 55
     return mQueue[n];
56 56
 }
... ...
@@ -142,24 +142,24 @@ bool ProposalQueue::birth(AtomicDomain &domain)
142 142
 
143 143
 bool ProposalQueue::death(AtomicDomain &domain)
144 144
 {
145
-    Atom a = domain.randomAtom();
146
-    if (mUsedIndices.count(a.pos / mDimensionSize))
145
+    Atom* a = domain.randomAtom();
146
+    if (mUsedIndices.count(a->pos / mDimensionSize))
147 147
     {
148 148
         return false; // matrix conflict - can't compute gibbs mass or deltaLL
149 149
     }
150 150
 
151
-    mQueue.push_back(AtomicProposal('D', a.pos, a.mass));
152
-    mUsedIndices.insert(a.pos / mDimensionSize);
153
-    mUsedPositions.insert(a.pos);
151
+    mQueue.push_back(AtomicProposal('D', a));
152
+    mUsedIndices.insert(a->pos / mDimensionSize);
153
+    mUsedPositions.insert(a->pos);
154 154
     --mMinAtoms;
155 155
     return true;
156 156
 }
157 157
 
158 158
 bool ProposalQueue::move(AtomicDomain &domain)
159 159
 {
160
-    Atom a = domain.randomAtom();
161
-    uint64_t lbound = domain.hasLeft(a) ? domain.left(a).pos : 0;
162
-    uint64_t rbound = domain.hasRight(a) ? domain.right(a).pos : mDomainSize;
160
+    AtomNeighborhood hood = domain.randomAtomWithNeighbors();
161
+    uint64_t lbound = hood.hasLeft() ? hood.left->pos : 1;
162
+    uint64_t rbound = hood.hasRight() ? hood.right->pos : mDomainSize;
163 163
 
164 164
     if (!mUsedPositions.isEmptyInterval(lbound, rbound))
165 165
     {
... ...
@@ -167,54 +167,55 @@ bool ProposalQueue::move(AtomicDomain &domain)
167 167
     }
168 168
 
169 169
     uint64_t newLocation = mRng.uniform64(lbound, rbound - 1);
170
-    if (mUsedIndices.count(a.pos / mDimensionSize) || mUsedIndices.count(newLocation / mDimensionSize))
170
+    if (mUsedIndices.count(hood.center->pos / mDimensionSize) || mUsedIndices.count(newLocation / mDimensionSize))
171 171
     {
172 172
         return false; // matrix conflict - can't compute deltaLL
173 173
     }
174 174
 
175
-    mQueue.push_back(AtomicProposal('M', a.pos, a.mass, newLocation));
176
-    mUsedIndices.insert(a.pos / mDimensionSize);
175
+    mQueue.push_back(AtomicProposal('M', hood.center, newLocation));
176
+    mUsedIndices.insert(hood.center->pos / mDimensionSize);
177 177
     mUsedIndices.insert(newLocation / mDimensionSize);
178
-    mUsedPositions.insert(a.pos);
178
+    mUsedPositions.insert(hood.center->pos);
179 179
     mUsedPositions.insert(newLocation);
180 180
     return true;
181 181
 }
182 182
 
183 183
 bool ProposalQueue::exchange(AtomicDomain &domain)
184 184
 {
185
-    Atom a1 = domain.randomAtom();
186
-    Atom a2 = domain.hasRight(a1) ? domain.right(a1) : domain.front();
185
+    AtomNeighborhood hood = domain.randomAtomWithRightNeighbor();
186
+    Atom* a1 = hood.center;    
187
+    Atom* a2 = hood.hasRight() ? hood.right : domain.front();
187 188
 
188
-    if (domain.hasRight(a1)) // has neighbor
189
+    if (hood.hasRight()) // has neighbor
189 190
     {
190
-        if (!mUsedPositions.isEmptyInterval(a1.pos, a2.pos))
191
+        if (!mUsedPositions.isEmptyInterval(a1->pos, a2->pos))
191 192
         {
192 193
             return false;
193 194
         }
194 195
     }
195 196
     else // exchange with first atom
196 197
     {
197
-        if (!mUsedPositions.isEmptyInterval(a1.pos, mDomainSize))
198
+        if (!mUsedPositions.isEmptyInterval(a1->pos, mDomainSize))
198 199
         {
199 200
             return false;
200 201
         }
201 202
         
202
-        if (!mUsedPositions.isEmptyInterval(0, domain.front().pos))
203
+        if (!mUsedPositions.isEmptyInterval(0, domain.front()->pos))
203 204
         {
204 205
             return false;
205 206
         }
206 207
     }
207 208
 
208
-    if (mUsedIndices.count(a1.pos / mDimensionSize) || mUsedIndices.count(a2.pos / mDimensionSize))
209
+    if (mUsedIndices.count(a1->pos / mDimensionSize) || mUsedIndices.count(a2->pos / mDimensionSize))
209 210
     {
210 211
         return false; // matrix conflict - can't compute gibbs mass or deltaLL
211 212
     }
212 213
 
213
-    mQueue.push_back(AtomicProposal('E', a1.pos, a1.mass, a2.pos, a2.mass));
214
-    mUsedIndices.insert(a1.pos / mDimensionSize);
215
-    mUsedIndices.insert(a2.pos / mDimensionSize);
216
-    mUsedPositions.insert(a1.pos);
217
-    mUsedPositions.insert(a2.pos);
214
+    mQueue.push_back(AtomicProposal('E', a1, a2));
215
+    mUsedIndices.insert(a1->pos / mDimensionSize);
216
+    mUsedIndices.insert(a2->pos / mDimensionSize);
217
+    mUsedPositions.insert(a1->pos);
218
+    mUsedPositions.insert(a2->pos);
218 219
     --mMinAtoms;
219 220
     return true;
220 221
 }
... ...
@@ -14,14 +14,32 @@
14 14
 struct AtomicProposal
15 15
 {
16 16
     char type;
17
-    uint64_t pos1;
18
-    float mass1;
19
-    uint64_t pos2;
20
-    float mass2;
17
+    uint64_t birthPos; // used in birth
18
+    uint64_t moveDest; // used in move
19
+
20
+    Atom *atom1; // used in death, move, exchange
21
+    Atom *atom2; // used in exchange
22
+
21 23
     mutable GapsRng rng;
22 24
 
23
-    AtomicProposal(char t, uint64_t p1, float m1=0.f, uint64_t p2=0, float m2=0.f)
24
-        : type(t), pos1(p1), mass1(m1), pos2(p2), mass2(m2)
25
+    // birth
26
+    AtomicProposal(char t, uint64_t pos)
27
+        : type(t), birthPos(pos), moveDest(0), atom1(NULL), atom2(NULL)
28
+    {}
29
+        
30
+    // death
31
+    AtomicProposal(char t, Atom *atom)
32
+        : type(t), birthPos(0), moveDest(0), atom1(atom), atom2(NULL)
33
+    {}
34
+
35
+    // move
36
+    AtomicProposal(char t, Atom *atom, uint64_t dest)
37
+        : type(t), birthPos(0), moveDest(dest), atom1(atom), atom2(NULL)
38
+    {}
39
+
40
+    // exchange
41
+    AtomicProposal(char t, Atom *a1, Atom *a2)
42
+        : type(t), birthPos(0), moveDest(0), atom1(a1), atom2(a2)
25 43
     {}
26 44
 };
27 45
 
... ...
@@ -74,7 +92,7 @@ public:
74 92
     void populate(AtomicDomain &domain, unsigned limit);
75 93
     void clear();
76 94
     unsigned size() const;
77
-    const AtomicProposal& operator[](int n) const;
95
+    AtomicProposal& operator[](int n);
78 96
 
79 97
     // update min/max atoms
80 98
     void acceptDeath();