Browse code

use a sorted vector for atomic domain

Tom Sherman authored on 23/08/2018 17:53:21
Showing14 changed files

... ...
@@ -1,5 +1,5 @@
1 1
 Package: CoGAPS
2
-Version: 3.3.11
2
+Version: 3.3.21
3 3
 Date: 2018-04-24
4 4
 Title: Coordinated Gene Activity in Pattern Sets
5 5
 Author: Thomas Sherman, Wai-shing Lee, Conor Kelton, Ondrej Maxian, Jacob Carey,
... ...
@@ -1,4 +1,4 @@
1
-# CoGAPS Version: 3.3.11
1
+# CoGAPS Version: 3.3.21
2 2
 
3 3
 [![Bioc](https://bioconductor.org/images/logo_bioconductor.gif)](https://bioconductor.org/packages/CoGAPS)
4 4
 [![downloads](https://bioconductor.org/shields/downloads/CancerInSilico.svg)](https://bioconductor.org/packages/CoGAPS)
... ...
@@ -1,6 +1,6 @@
1 1
 #! /bin/sh
2 2
 # Guess values for system-dependent variables and create Makefiles.
3
-# Generated by GNU Autoconf 2.69 for CoGAPS 3.3.2.
3
+# Generated by GNU Autoconf 2.69 for CoGAPS 3.3.11.
4 4
 #
5 5
 #
6 6
 # Copyright (C) 1992-1996, 1998-2012 Free Software Foundation, Inc.
... ...
@@ -577,8 +577,8 @@ MAKEFLAGS=
577 577
 # Identity of this package.
578 578
 PACKAGE_NAME='CoGAPS'
579 579
 PACKAGE_TARNAME='cogaps'
580
-PACKAGE_VERSION='3.3.2'
581
-PACKAGE_STRING='CoGAPS 3.3.2'
580
+PACKAGE_VERSION='3.3.11'
581
+PACKAGE_STRING='CoGAPS 3.3.11'
582 582
 PACKAGE_BUGREPORT=''
583 583
 PACKAGE_URL=''
584 584
 
... ...
@@ -696,6 +696,7 @@ enable_option_checking
696 696
 enable_silent_rules
697 697
 enable_dependency_tracking
698 698
 enable_debug
699
+enable_internal_tests
699 700
 enable_warnings
700 701
 enable_simd
701 702
 enable_openmp
... ...
@@ -1262,7 +1263,7 @@ if test "$ac_init_help" = "long"; then
1262 1263
   # Omit some internal or obsolete options to make the list less imposing.
1263 1264
   # This message is too long to be a string in the A/UX 3.1 sh.
1264 1265
   cat <<_ACEOF
1265
-\`configure' configures CoGAPS 3.3.2 to adapt to many kinds of systems.
1266
+\`configure' configures CoGAPS 3.3.11 to adapt to many kinds of systems.
1266 1267
 
1267 1268
 Usage: $0 [OPTION]... [VAR=VALUE]...
1268 1269
 
... ...
@@ -1333,7 +1334,7 @@ fi
1333 1334
 
1334 1335
 if test -n "$ac_init_help"; then
1335 1336
   case $ac_init_help in
1336
-     short | recursive ) echo "Configuration of CoGAPS 3.3.2:";;
1337
+     short | recursive ) echo "Configuration of CoGAPS 3.3.11:";;
1337 1338
    esac
1338 1339
   cat <<\_ACEOF
1339 1340
 
... ...
@@ -1348,6 +1349,7 @@ Optional Features:
1348 1349
   --disable-dependency-tracking
1349 1350
                           speeds up one-time build
1350 1351
   --enable-debug          build debug version of CoGAPS
1352
+  --enable-internal-tests allow for testing private data
1351 1353
   --enable-warnings       compile CoGAPS with warning messages
1352 1354
   --enable-simd           specify simd instruction set (sse, avx)
1353 1355
   --enable-openmp         compile with openMP support if available
... ...
@@ -1430,7 +1432,7 @@ fi
1430 1432
 test -n "$ac_init_help" && exit $ac_status
1431 1433
 if $ac_init_version; then
1432 1434
   cat <<\_ACEOF
1433
-CoGAPS configure 3.3.2
1435
+CoGAPS configure 3.3.11
1434 1436
 generated by GNU Autoconf 2.69
1435 1437
 
1436 1438
 Copyright (C) 2012 Free Software Foundation, Inc.
... ...
@@ -1873,7 +1875,7 @@ cat >config.log <<_ACEOF
1873 1875
 This file contains any messages produced by compilers while
1874 1876
 running configure, to aid debugging if configure makes a mistake.
1875 1877
 
1876
-It was created by CoGAPS $as_me 3.3.2, which was
1878
+It was created by CoGAPS $as_me 3.3.11, which was
1877 1879
 generated by GNU Autoconf 2.69.  Invocation command line was
1878 1880
 
1879 1881
   $ $0 $@
... ...
@@ -2736,7 +2738,7 @@ fi
2736 2738
 
2737 2739
 # Define the identity of the package.
2738 2740
  PACKAGE='cogaps'
2739
- VERSION='3.3.2'
2741
+ VERSION='3.3.11'
2740 2742
 
2741 2743
 
2742 2744
 cat >>confdefs.h <<_ACEOF
... ...
@@ -4069,6 +4071,15 @@ else
4069 4071
 fi
4070 4072
 
4071 4073
 
4074
+# Check if running internal tests
4075
+# Check whether --enable-internal-tests was given.
4076
+if test "${enable_internal_tests+set}" = set; then :
4077
+  enableval=$enable_internal_tests; internal_tests=yes
4078
+else
4079
+  internal_tests=no
4080
+fi
4081
+
4082
+
4072 4083
 # Check if compiler warnings should be turned on
4073 4084
 # Check whether --enable-warnings was given.
4074 4085
 if test "${enable_warnings+set}" = set; then :
... ...
@@ -5140,7 +5151,6 @@ else
5140 5151
 fi
5141 5152
 
5142 5153
 
5143
-if test "x$use_openmp" != "xno" ; then
5144 5154
 
5145 5155
 
5146 5156
 { $as_echo "$as_me:${as_lineno-$LINENO}: checking for OpenMP flag of C++ compiler" >&5
... ...
@@ -5207,6 +5217,7 @@ $as_echo "#define HAVE_OPENMP 1" >>confdefs.h
5207 5217
 
5208 5218
 fi
5209 5219
 
5220
+if test "x$use_openmp" != "xno" ; then
5210 5221
     GAPS_CXX_FLAGS+=" $OPENMP_CXXFLAGS "
5211 5222
     GAPS_LIBS+=" $OPENMP_CXXFLAGS "
5212 5223
 fi
... ...
@@ -6841,6 +6852,12 @@ if test "x$build_debug" = "xyes" ; then
6841 6852
     GAPS_CPP_FLAGS+=" -DGAPS_DEBUG "
6842 6853
 fi
6843 6854
 
6855
+# set compile flags for internal tests
6856
+if test "x$internal_tests" = "xyes" ; then
6857
+    echo "Enabling Internal Tests"
6858
+    GAPS_CPP_FLAGS+=" -DGAPS_INTERNAL_TESTS "
6859
+fi
6860
+
6844 6861
 # set compile flags if warnings enabled
6845 6862
 if test "x$warnings" = "xyes" ; then
6846 6863
     if test "x$ax_cv_cxx_compiler_vendor" = "xgnu" || test "x$ax_cv_cxx_compiler_vendor" = "xclang"; then
... ...
@@ -7472,7 +7489,7 @@ cat >>$CONFIG_STATUS <<\_ACEOF || ac_write_fail=1
7472 7489
 # report actual input values of CONFIG_FILES etc. instead of their
7473 7490
 # values after options handling.
7474 7491
 ac_log="
7475
-This file was extended by CoGAPS $as_me 3.3.2, which was
7492
+This file was extended by CoGAPS $as_me 3.3.11, which was
7476 7493
 generated by GNU Autoconf 2.69.  Invocation command line was
7477 7494
 
7478 7495
   CONFIG_FILES    = $CONFIG_FILES
... ...
@@ -7529,7 +7546,7 @@ _ACEOF
7529 7546
 cat >>$CONFIG_STATUS <<_ACEOF || ac_write_fail=1
7530 7547
 ac_cs_config="`$as_echo "$ac_configure_args" | sed 's/^ //; s/[\\""\`\$]/\\\\&/g'`"
7531 7548
 ac_cs_version="\\
7532
-CoGAPS config.status 3.3.2
7549
+CoGAPS config.status 3.3.11
7533 7550
 configured by $0, generated by GNU Autoconf 2.69,
7534 7551
   with options \\"\$ac_cs_config\\"
7535 7552
 
... ...
@@ -14,6 +14,10 @@ AC_PROG_CXX
14 14
 AC_ARG_ENABLE(debug, [AC_HELP_STRING([--enable-debug],
15 15
     [build debug version of CoGAPS])], [build_debug=yes], [build_debug=no])
16 16
 
17
+# Check if running internal tests
18
+AC_ARG_ENABLE(internal-tests, [AC_HELP_STRING([--enable-internal-tests],
19
+    [allow for testing private data])], [internal_tests=yes], [internal_tests=no])
20
+
17 21
 # Check if compiler warnings should be turned on
18 22
 AC_ARG_ENABLE(warnings, [AC_HELP_STRING([--enable-warnings],
19 23
     [compile CoGAPS with warning messages])], [warnings=yes], [warnings=no])
... ...
@@ -54,6 +58,12 @@ if test "x$build_debug" = "xyes" ; then
54 58
     GAPS_CPP_FLAGS+=" -DGAPS_DEBUG "
55 59
 fi
56 60
 
61
+# set compile flags for internal tests
62
+if test "x$internal_tests" = "xyes" ; then
63
+    echo "Enabling Internal Tests"
64
+    GAPS_CPP_FLAGS+=" -DGAPS_INTERNAL_TESTS "
65
+fi
66
+
57 67
 # set compile flags if warnings enabled
58 68
 if test "x$warnings" = "xyes" ; then
59 69
     if test "x$ax_cv_cxx_compiler_vendor" = "xgnu" || test "x$ax_cv_cxx_compiler_vendor" = "xclang"; then
... ...
@@ -2,8 +2,7 @@
2 2
 #define __COGAPS_ARCHIVE_H__
3 3
 
4 4
 #include <fstream>
5
-
6
-#include <boost/random/mersenne_twister.hpp>
5
+#include <stdint.h>
7 6
 
8 7
 // flags for opening an archive
9 8
 #define ARCHIVE_READ  std::ios::in
... ...
@@ -53,7 +52,6 @@ public:
53 52
     friend Archive& operator<<(Archive &ar, int64_t val)  { return writeToArchive(ar, val); }
54 53
     friend Archive& operator<<(Archive &ar, float val)    { return writeToArchive(ar, val); }
55 54
     friend Archive& operator<<(Archive &ar, double val)   { return writeToArchive(ar, val); }
56
-    friend Archive& operator<<(Archive &ar, boost::random::mt11213b val) { return writeToArchive(ar, val); }
57 55
 
58 56
     friend Archive& operator>>(Archive &ar, char &val)     { return readFromArchive(ar, val); }
59 57
     friend Archive& operator>>(Archive &ar, bool &val)     { return readFromArchive(ar, val); }
... ...
@@ -63,7 +61,6 @@ public:
63 61
     friend Archive& operator>>(Archive &ar, int64_t &val)  { return readFromArchive(ar, val); }
64 62
     friend Archive& operator>>(Archive &ar, float &val)    { return readFromArchive(ar, val); }
65 63
     friend Archive& operator>>(Archive &ar, double &val)   { return readFromArchive(ar, val); }
66
-    friend Archive& operator>>(Archive &ar, boost::random::mt11213b &val) { return readFromArchive(ar, val); }
67 64
 };
68 65
 
69 66
 #endif // __COGAPS_ARCHIVE_H__
... ...
@@ -3,87 +3,6 @@
3 3
 
4 4
 #include <vector>
5 5
 
6
-// should only ever use this when copying a default constructed bucket
7
-// in the vector of buckets used in the allocator
8
-/*AtomBucket& AtomBucket::operator=(const AtomBucket& other)
9
-{
10
-    GAPS_ASSERT(other.mSize == 0);
11
-    GAPS_ASSERT(other.mPrev == NULL);
12
-    GAPS_ASSERT(other.mNext == NULL);
13
-    GAPS_ASSERT(other.mOverflow == NULL);
14
-
15
-    mSize = 0;
16
-    mPrev = NULL;
17
-    mNext = NULL;
18
-    mOverflow = NULL;
19
-}*/
20
-
21
-/////////////////////// ALLOCATOR FOR OVERFLOW BUCKETS /////////////////////////
22
-
23
-static const unsigned poolsize = 1024;
24
-
25
-class AtomBucketPool
26
-{
27
-public:
28
-
29
-    AtomBucketPool()
30
-        : mPool(std::vector<AtomBucket>(poolsize, AtomBucket())), mCurrent(0)
31
-    {}
32
-
33
-    AtomBucket* create()
34
-    {
35
-        return &(mPool[mCurrent++]);
36
-    }
37
-
38
-    bool depleted()
39
-    {
40
-        return mCurrent == poolsize;
41
-    }
42
-
43
-private:
44
-
45
-    std::vector<AtomBucket> mPool;
46
-    unsigned mCurrent;
47
-};
48
-
49
-class AtomBucketAllocator
50
-{
51
-public:
52
-
53
-    AtomBucketAllocator()
54
-    {
55
-        mAllPools.push_back(new AtomBucketPool());
56
-    }
57
-
58
-    ~AtomBucketAllocator()
59
-    {
60
-        for (unsigned i = 0; i < mAllPools.size(); ++i)
61
-        {
62
-            delete mAllPools[i];
63
-        }
64
-    }
65
-
66
-    AtomBucket* create()
67
-    {
68
-        if (mAllPools.back()->depleted())
69
-        {
70
-            mAllPools.push_back(new AtomBucketPool());
71
-        }
72
-        return mAllPools.back()->create();
73
-    }
74
-
75
-private:
76
-
77
-    std::vector<AtomBucketPool*> mAllPools;
78
-};
79
-
80
-static AtomBucket* createAtomBucket()
81
-{
82
-    //static AtomBucketAllocator allocator;
83
-    //return allocator.create();
84
-    return new AtomBucket();
85
-}
86
-
87 6
 ////////////////////////////////// ATOM ////////////////////////////////////////
88 7
 
89 8
 Atom::Atom() : pos(0), mass(0.f) {}
... ...
@@ -125,583 +44,71 @@ bool AtomNeighborhood::hasRight()
125 44
     return right != NULL;
126 45
 }
127 46
 
128
-/////////////////////////////// ATOM BUCKET ////////////////////////////////////
129
-
130
-// note much of this code is unrolled intentionally so that the most frequent
131
-// cases are fast
132
-
133
-AtomBucket::AtomBucket()
134
-    : mSize(0), mOverflow(NULL), mPrev(NULL), mNext(NULL)
135
-{}
136
-
137
-unsigned AtomBucket::size() const
138
-{
139
-    GAPS_ASSERT(mPrev == NULL || !mPrev->isEmpty());
140
-    GAPS_ASSERT(mNext == NULL || !mNext->isEmpty());
141
-
142
-    return mOverflow == NULL ? mSize : mSize + mOverflow->size();
143
-}
144
-
145
-bool AtomBucket::isEmpty() const
146
-{
147
-    return mSize == 0;
148
-}
149
-
150
-bool AtomBucket::contains(uint64_t pos) const
151
-{
152
-    GAPS_ASSERT(mPrev == NULL || !mPrev->isEmpty());
153
-    GAPS_ASSERT(mNext == NULL || !mNext->isEmpty());
154
-    GAPS_ASSERT(pos > 0);
155
-
156
-    if (mOverflow != NULL && !mOverflow->isEmpty() && pos > mBucket[1].pos)
157
-    {
158
-        return mOverflow->contains(pos);
159
-    }
160
-    else
161
-    {
162
-        switch (mSize)
163
-        {
164
-            case 0:
165
-                return false;
166
-            case 1:
167
-                return mBucket[0].pos == pos;
168
-            case 2:
169
-                return mBucket[0].pos == pos || mBucket[1].pos == pos;
170
-        }
171
-    }
172
-}
173
-
174
-Atom* AtomBucket::operator[](unsigned index)
175
-{
176
-    GAPS_ASSERT(mPrev == NULL || !mPrev->isEmpty());
177
-    GAPS_ASSERT(mNext == NULL || !mNext->isEmpty());
178
-    GAPS_ASSERT(index < size());
179
-
180
-    if (index > 1) // must be overflowed
181
-    {
182
-        return mOverflow->operator[](index - 2);
183
-    }
184
-    return &(mBucket[index]);
185
-}
186
-
187
-void AtomBucket::insert(uint64_t pos, float mass)
188
-{
189
-    GAPS_ASSERT(mPrev == NULL || !mPrev->isEmpty());
190
-    GAPS_ASSERT(mNext == NULL || !mNext->isEmpty());
191
-    GAPS_ASSERT(pos > 0);
192
-
193
-    if (mSize == 0)
194
-    {
195
-        mBucket[0] = Atom(pos, mass);
196
-        ++mSize;
197
-    }
198
-    else if (mSize == 1)
199
-    {
200
-        if (pos < mBucket[0].pos)
201
-        {
202
-            mBucket[1] = mBucket[0];
203
-            mBucket[0] = Atom(pos, mass);
204
-        }
205
-        else
206
-        {
207
-            mBucket[1] = Atom(pos, mass);
208
-        }
209
-        ++mSize;
210
-    }
211
-    else
212
-    {
213
-        // check if we need to allocate the overflow bucket
214
-        if (mOverflow == NULL)
215
-        {
216
-            mOverflow = createAtomBucket();
217
-            mOverflow->mPrev = this;
218
-            mOverflow->mNext = mNext;
219
-        }
220
-        else if (mOverflow->isEmpty())
221
-        {
222
-            mOverflow->mPrev = this;
223
-            mOverflow->mNext = mNext;
224
-        }
225
-
226
-        // push correct atom into overflow bucket
227
-        if (pos > mBucket[1].pos)
228
-        {
229
-            return mOverflow->insert(pos, mass);
230
-        }
231
-        mOverflow->insert(mBucket[1].pos, mBucket[1].mass);
232
-
233
-        // if inserting in this bucket, find correct position
234
-        if (pos < mBucket[0].pos)
235
-        {
236
-            mBucket[1] = mBucket[0];
237
-            mBucket[0] = Atom(pos, mass);
238
-        }
239
-        else
240
-        {
241
-            mBucket[1] = Atom(pos, mass);
242
-        }
243
-    }
244
-}
245
-
246
-// assumes pos is contained in this chain
247
-void AtomBucket::erase(uint64_t pos)
248
-{
249
-    GAPS_ASSERT(mPrev == NULL || !mPrev->isEmpty());
250
-    GAPS_ASSERT(mNext == NULL || !mNext->isEmpty());
251
-    GAPS_ASSERT(pos > 0);
252
-    GAPS_ASSERT(mSize > 0);
253
-    GAPS_ASSERT(contains(pos));
254
-
255
-    if (mSize == 1)
256
-    {
257
-        connectAdjacent();
258
-        mSize = 0;
259
-    }
260
-    else if (mSize == 2)
261
-    {
262
-        // check if this position is in overflow bucket    
263
-        if (pos > mBucket[1].pos)
264
-        {
265
-            return mOverflow->erase(pos);
266
-        }
267
-
268
-        // shift top position down if needed
269
-        if (mBucket[0].pos == pos)
270
-        {
271
-            mBucket[0] = mBucket[1];
272
-        }
273
-        
274
-        // pull first atom from overflow if it's there
275
-        if (mOverflow != NULL && !mOverflow->isEmpty())
276
-        {
277
-            mBucket[1] = *(mOverflow->front());
278
-            mOverflow->eraseFront();
279
-        }
280
-        else // just delete atom at last position
281
-        {
282
-            --mSize;
283
-        }
284
-    }
285
-}
286
-
287
-void AtomBucket::eraseFront()
288
-{
289
-    GAPS_ASSERT(mPrev == NULL || !mPrev->isEmpty());
290
-    GAPS_ASSERT(mNext == NULL || !mNext->isEmpty());
291
-
292
-    GAPS_ASSERT(mSize > 0);
293
-
294
-    if (mSize == 1)
295
-    {
296
-        connectAdjacent();
297
-        mSize = 0;
298
-    }
299
-    else if (mSize == 2)
300
-    {
301
-        mBucket[0] = mBucket[1];
302
-        
303
-        // pull first atom from overflow if it's there
304
-        if (mOverflow != NULL && !mOverflow->isEmpty())
305
-        {
306
-            mBucket[1] = *(mOverflow->front());
307
-            mOverflow->eraseFront();
308
-        }
309
-        else
310
-        {
311
-            --mSize;
312
-        }
313
-    }
314
-}
315
-
316
-AtomNeighborhood AtomBucket::getNeighbors(unsigned index)
317
-{
318
-    GAPS_ASSERT(mPrev == NULL || !mPrev->isEmpty());
319
-    GAPS_ASSERT(mNext == NULL || !mNext->isEmpty());
320
-    GAPS_ASSERT(index < size());
321
-
322
-    return AtomNeighborhood(getLeft(index), this->operator[](index), getRight(index));
323
-}
324
-
325
-AtomNeighborhood AtomBucket::getRightNeighbor(unsigned index)
326
-{
327
-    GAPS_ASSERT(mPrev == NULL || !mPrev->isEmpty());
328
-    GAPS_ASSERT(mNext == NULL || !mNext->isEmpty());
329
-    GAPS_ASSERT(index < size());
330
-
331
-    return AtomNeighborhood(NULL, this->operator[](index), getRight(index));
332
-}
333
-
334
-// needs to propogate through overflow so that the last overflow will know
335
-// where the next bucket is
336
-void AtomBucket::setRightAdjacentBucket(AtomBucket *bucket)
337
-{
338
-    GAPS_ASSERT(mPrev == NULL || !mPrev->isEmpty());
339
-    GAPS_ASSERT(mNext == NULL || !mNext->isEmpty());
340
-
341
-    mNext = bucket;
342
-    if (mOverflow != NULL && !mOverflow->isEmpty())
343
-    {
344
-        mOverflow->setRightAdjacentBucket(bucket);
345
-    }
346
-}
347
-
348
-void AtomBucket::setLeftAdjacentBucket(AtomBucket *bucket)
349
-{
350
-    GAPS_ASSERT(mPrev == NULL || !mPrev->isEmpty());
351
-    GAPS_ASSERT(mNext == NULL || !mNext->isEmpty());
352
-    GAPS_ASSERT(bucket == NULL || !bucket->isEmpty());
353
-
354
-    mPrev = bucket;
355
-}
356
-
357
-void AtomBucket::connectAdjacent()
358
-{
359
-    GAPS_ASSERT(mPrev == NULL || !mPrev->isEmpty());
360
-    GAPS_ASSERT(mNext == NULL || !mNext->isEmpty());
361
-
362
-    // be careful when connecting overflow buckets
363
-    if (mNext != NULL && (mPrev == NULL || mPrev->mOverflow != this))
364
-    {
365
-        mNext->setLeftAdjacentBucket(mPrev);
366
-    }
367
-    if (mPrev != NULL)
368
-    {
369
-        mPrev->setRightAdjacentBucket(mNext);
370
-    }
371
-
372
-    mNext = NULL;
373
-    mPrev = NULL;
374
-}
375
-
376
-Atom* AtomBucket::front()
377
-{
378
-    GAPS_ASSERT(mPrev == NULL || !mPrev->isEmpty());
379
-    GAPS_ASSERT(mNext == NULL || !mNext->isEmpty());
380
-    GAPS_ASSERT(mSize > 0);
381
-    GAPS_ASSERT(mBucket[0].pos > 0);
382
-
383
-    return &(mBucket[0]);
384
-}
385
-
386
-Atom* AtomBucket::back()
387
-{
388
-    GAPS_ASSERT(mPrev == NULL || !mPrev->isEmpty());
389
-    GAPS_ASSERT(mNext == NULL || !mNext->isEmpty());
390
-    GAPS_ASSERT(mSize > 0);
391
-
392
-    if (mOverflow == NULL || (mOverflow != NULL && mOverflow->isEmpty()))
393
-    {
394
-        return &(mBucket[mSize - 1]);
395
-    }
396
-    return mOverflow->back();
397
-}
398
-
399
-Atom* AtomBucket::getLeft(unsigned index)
400
-{
401
-    GAPS_ASSERT(mPrev == NULL || !mPrev->isEmpty());
402
-    GAPS_ASSERT(mNext == NULL || !mNext->isEmpty());
403
-    GAPS_ASSERT(mSize > 0);
404
-    GAPS_ASSERT(index < size());
405
-
406
-    if (index == 0)
407
-    {
408
-        GAPS_ASSERT(mPrev == NULL || mPrev->mOverflow != this);
409
-        return mPrev != NULL ? mPrev->back() : NULL; // mPrev can't be overflow here
410
-    }
411
-    else if (index < 3)
412
-    {
413
-        return &(mBucket[index - 1]);
414
-    }
415
-    else
416
-    {
417
-        GAPS_ASSERT(mOverflow != NULL);
418
-        return mOverflow->getLeft(index - 2);
419
-    }    
420
-}
421
-
422
-Atom* AtomBucket::getRight(unsigned index)
423
-{
424
-    GAPS_ASSERT(mPrev == NULL || !mPrev->isEmpty());
425
-    GAPS_ASSERT(mNext == NULL || !mNext->isEmpty());
426
-    GAPS_ASSERT(mSize > 0);
427
-    GAPS_ASSERT(index < size());
428
-
429
-    if (index == mSize - 1) // this is last atom in bucket
430
-    {
431
-        if (mOverflow != NULL && !mOverflow->isEmpty())
432
-        {
433
-            return mOverflow->front();
434
-        }
435
-        return mNext != NULL ? mNext->front() : NULL;
436
-    }        
437
-    else if (index == 0)
438
-    {
439
-        GAPS_ASSERT(mBucket[1].pos > 0);
440
-        return &(mBucket[1]);
441
-    }
442
-    else
443
-    {
444
-        mOverflow->getRight(index - 2);
445
-    }
446
-}
447
-
448
-Archive& operator<<(Archive &ar, AtomBucket &b)
449
-{
450
-    bool hasOverflow = (b.mOverflow != NULL);
451
-    ar << b.mBucket[0] << b.mBucket[1] << b.mSize << hasOverflow;
452
-
453
-    if (hasOverflow)
454
-    {
455
-        ar << *b.mOverflow;
456
-    }
457
-    return ar;
458
-}
459
-
460
-Archive& operator>>(Archive& ar, AtomBucket &b)
461
-{
462
-    bool hasOverflow = false;
463
-    ar >> b.mBucket[0] >> b.mBucket[1] >> b.mSize >> hasOverflow;
464
-
465
-    if (hasOverflow)
466
-    {
467
-        b.mOverflow = createAtomBucket();
468
-        ar >> *b.mOverflow;
469
-    }
470
-    return ar;
471
-}
472
-
473
-////////////////////////////// ATOM HASH MAP ///////////////////////////////////
474
-
475
-AtomHashMap::AtomHashMap(unsigned expectedNumAtoms)
476
-    :
477
-mExpectedNumAtoms(expectedNumAtoms), mLongestBucketSize(0), mSize(0),
478
-mHashMap(std::vector<AtomBucket>(expectedNumAtoms, AtomBucket()))
479
-{}
480
-
481
-void AtomHashMap::setTotalLength(uint64_t length)
482
-{
483
-    GAPS_ASSERT(length % mExpectedNumAtoms == 0);
484
-    mTotalLength = length;
485
-    mBucketLength = mTotalLength / mExpectedNumAtoms;
486
-}
487
-
488
-// pos ranges from 1 to mTotalLength
489
-unsigned AtomHashMap::hash(uint64_t pos) const
490
-{
491
-    return pos / mBucketLength;
492
-}
493
-
494
-Atom* AtomHashMap::front()
495
-{
496
-    return mHashMap[*(mFullBuckets.begin())].front();
497
-}
498
-    
499
-HashMapIndex AtomHashMap::getRandomIndex() const
500
-{
501
-    GAPS_ASSERT(mSize > 0);
502
-    while (true)
503
-    {
504
-        unsigned bucket = hash(mRng.uniform64(0, mTotalLength - 1));
505
-        if (!mHashMap[bucket].isEmpty())
506
-        {
507
-            unsigned pos = mRng.uniform32(0, mLongestBucketSize - 1);
508
-            if (pos < mHashMap[bucket].size())
509
-            {
510
-                return HashMapIndex(bucket, pos);
511
-            }
512
-        }
513
-    }
514
-}
515
-
516
-Atom* AtomHashMap::randomAtom()
517
-{
518
-    HashMapIndex ndx = getRandomIndex();
519
-    return mHashMap[ndx.bucket][ndx.index];
520
-}
521
-
522
-AtomNeighborhood AtomHashMap::randomAtomWithNeighbors()
523
-{
524
-    HashMapIndex ndx = getRandomIndex();
525
-    return mHashMap[ndx.bucket].getNeighbors(ndx.index);
526
-}
527
-
528
-AtomNeighborhood AtomHashMap::randomAtomWithRightNeighbor()
529
-{
530
-    HashMapIndex ndx = getRandomIndex();
531
-    return mHashMap[ndx.bucket].getRightNeighbor(ndx.index);
532
-}
533
-
534
-bool AtomHashMap::contains(uint64_t pos) const
535
-{
536
-    unsigned index = hash(pos);
537
-    return !(mHashMap[index].isEmpty()) && mHashMap[index].contains(pos);
538
-}
539
-
540
-unsigned AtomHashMap::size() const
541
-{
542
-    return mSize;
543
-}
544
-
545
-// usually O(1), might be O(logN)
546
-void AtomHashMap::insert(uint64_t pos, float mass)
547
-{
548
-    // hash position once
549
-    unsigned index = hash(pos);
550
-
551
-    // if inserting into an empty bucket, mark bucket as non-empty
552
-    // and connect adjacent buckets O(logN)
553
-    if (mHashMap[index].isEmpty())
554
-    {
555
-        mHashMap[index].insert(pos, mass);
556
-        GAPS_ASSERT(mHashMap[index].contains(pos));
557
-        std::set<unsigned>::iterator it = mFullBuckets.insert(index).first;
558
-        if (it != mFullBuckets.begin())
559
-        {
560
-            --it;
561
-            mHashMap[index].setLeftAdjacentBucket(&mHashMap[*it]);
562
-            mHashMap[*it].setRightAdjacentBucket(&mHashMap[index]);
563
-            ++it;
564
-        }
565
-        if (++it != mFullBuckets.end())
566
-        {
567
-            mHashMap[index].setRightAdjacentBucket(&mHashMap[*it]);
568
-            mHashMap[*it].setLeftAdjacentBucket(&mHashMap[index]);
569
-        }
570
-    }
571
-    else
572
-    {
573
-        mHashMap[index].insert(pos, mass);
574
-        GAPS_ASSERT(mHashMap[index].contains(pos));
575
-    }
576
-    GAPS_ASSERT(mHashMap[index].contains(pos));
577
-
578
-    // insert atom
579
-    ++mSize;
580
-
581
-    // check if this is now the longest bucket
582
-    if (mHashMap[index].size() > mLongestBucketSize)
583
-    {
584
-        mWhichLongestBucket = index;
585
-        ++mLongestBucketSize;
586
-    }
587
-}
588
-
589
-// usually O(1), sometimes O(logN), rarely O(N)
590
-void AtomHashMap::erase(uint64_t pos)
591
-{
592
-    // erase atom
593
-    unsigned index = hash(pos);
594
-    GAPS_ASSERT(!mHashMap[index].isEmpty());
595
-    GAPS_ASSERT(mHashMap[index].contains(pos));
596
-    --mSize;
597
-    
598
-    // mark bucket as empty if this was the last atom O(logN)
599
-    mHashMap[index].erase(pos);
600
-    if (mHashMap[index].isEmpty())
601
-    {
602
-        mFullBuckets.erase(index);
603
-    }
604
-    GAPS_ASSERT(!mHashMap[index].contains(pos));
605
-
606
-    // if this atom was in the largest bucket, find the new largest bucket O(N)
607
-    if (index == mWhichLongestBucket)
608
-    {
609
-        --mLongestBucketSize;
610
-        std::set<unsigned>::iterator it = mFullBuckets.begin();
611
-        for (; it != mFullBuckets.end(); ++it)
612
-        {
613
-            if (mHashMap[*it].size() > mLongestBucketSize)
614
-            {
615
-                mLongestBucketSize = mHashMap[*it].size();
616
-                mWhichLongestBucket = *it;
617
-            }
618
-        }
619
-    }
620
-}
47
+////////////////////////////// ATOMIC DOMAIN ///////////////////////////////////
621 48
 
622
-void AtomHashMap::move(uint64_t src, uint64_t dest, float mass)
49
+AtomicDomain::AtomicDomain(unsigned nBins)
623 50
 {
624
-    unsigned srcHash = hash(src);
625
-    unsigned destHash = hash(dest);
626
-    if (srcHash != destHash)
627
-    {
628
-        erase(src);
629
-        insert(dest, mass);
630
-    }
51
+    uint64_t binLength = std::numeric_limits<uint64_t>::max() / nBins;
52
+    mDomainLength = binLength * nBins;
631 53
 }
632 54
 
633
-Archive& operator<<(Archive& ar, AtomHashMap &h)
55
+Atom* AtomicDomain::front()
634 56
 {
635
-    ar << h.mExpectedNumAtoms << h.mBucketLength << h.mTotalLength <<
636
-        h.mWhichLongestBucket << h.mLongestBucketSize << h.mSize;
637
-
638
-    ar << h.mFullBuckets.size();
639
-    std::set<unsigned>::iterator it = h.mFullBuckets.begin();
640
-    for (; it != h.mFullBuckets.end(); ++it)
641
-    {
642
-        ar << *it << h.mHashMap[*it];
643
-    }
644
-    return ar;
57
+    GAPS_ASSERT(size() > 0);
58
+    return &(mAtoms.front());
645 59
 }
646 60
 
647
-Archive& operator>>(Archive& ar, AtomHashMap &h)
61
+Atom* AtomicDomain::randomAtom()
648 62
 {
649
-    ar >> h.mExpectedNumAtoms >> h.mBucketLength >> h.mTotalLength >>
650
-        h.mWhichLongestBucket >> h.mLongestBucketSize >> h.mSize;
651
-
652
-    unsigned nBuckets = 0;
653
-    ar >> nBuckets;
654
-
655
-    for (unsigned i = 0; i < nBuckets; ++i)
656
-    {
657
-        unsigned thisBucket = 0;
658
-        ar >> thisBucket;
659
-        h.mFullBuckets.insert(thisBucket);
660
-        ar >> h.mHashMap[thisBucket];
661
-    }
662
-    return ar;
63
+    GAPS_ASSERT(size() > 0);
64
+    unsigned index = mRng.uniform32(0, mAtoms.size() - 1);
65
+    return &(mAtoms[index]);
663 66
 }
664 67
 
665
-////////////////////////////// ATOMIC DOMAIN ///////////////////////////////////
666
-
667
-AtomicDomain::AtomicDomain(unsigned nBins) : mAtoms(nBins)
68
+AtomNeighborhood AtomicDomain::randomAtomWithNeighbors()
668 69
 {
669
-    uint64_t binLength = std::numeric_limits<uint64_t>::max() / nBins;
670
-    mDomainLength = binLength * nBins;
671
-    mAtoms.setTotalLength(mDomainLength);
70
+    GAPS_ASSERT(size() > 0);
71
+    unsigned index = mRng.uniform32(0, mAtoms.size() - 1);
72
+    Atom* left = (index == 0) ? NULL : &(mAtoms[index - 1]);
73
+    Atom* right = (index == mAtoms.size() - 1) ? NULL : &(mAtoms[index + 1]);
74
+    return AtomNeighborhood(left, &(mAtoms[index]), right);
672 75
 }
673 76
 
674
-void AtomicDomain::setDomainSize(uint64_t size)
77
+AtomNeighborhood AtomicDomain::randomAtomWithRightNeighbor()
675 78
 {
676
-    // do nothing
79
+    GAPS_ASSERT(size() > 0);
80
+    unsigned index = mRng.uniform32(0, mAtoms.size() - 1);
81
+    Atom* right = (index == mAtoms.size() - 1) ? NULL : &(mAtoms[index + 1]);
82
+    return AtomNeighborhood(NULL, &(mAtoms[index]), right);
677 83
 }
678 84
 
679
-Atom* AtomicDomain::front()
85
+static bool compareAtomLower(const Atom &lhs, uint64_t pos)
680 86
 {
681
-    return mAtoms.front();
87
+    return lhs.pos < pos;
682 88
 }
683 89
 
684
-Atom* AtomicDomain::randomAtom()
90
+static bool compareAtomUpper(uint64_t pos, const Atom &lhs)
685 91
 {
686
-    return mAtoms.randomAtom();
92
+    return lhs.pos < pos;
687 93
 }
688 94
 
689
-AtomNeighborhood AtomicDomain::randomAtomWithNeighbors()
95
+static bool compareAtom(const Atom &lhs, const Atom &rhs)
690 96
 {
691
-    return mAtoms.randomAtomWithNeighbors();
97
+    return lhs.pos < rhs.pos;
692 98
 }
693 99
 
694
-AtomNeighborhood AtomicDomain::randomAtomWithRightNeighbor()
100
+static bool vecContains(const std::vector<Atom> &vec, uint64_t pos)
695 101
 {
696
-    return mAtoms.randomAtomWithRightNeighbor();
102
+    Atom temp(pos, 0.f);
103
+    return std::binary_search(vec.begin(), vec.end(), temp, compareAtom);
697 104
 }
698 105
 
699 106
 uint64_t AtomicDomain::randomFreePosition() const
700 107
 {
701
-    uint64_t pos = mRng.uniform64(1, mDomainLength);
702
-    while (mAtoms.contains(pos))
108
+    uint64_t pos = mRng.uniform64(0, mDomainLength);
109
+    while (vecContains(mAtoms, pos))
703 110
     {
704
-        pos = mRng.uniform64(1, mDomainLength);
111
+        pos = mRng.uniform64(0, mDomainLength);
705 112
     } 
706 113
     return pos;
707 114
 }
... ...
@@ -731,54 +138,65 @@ void AtomicDomain::cacheErase(uint64_t pos) const
731 138
     mEraseCache[ndx] = pos;
732 139
 }
733 140
 
734
-void AtomicDomain::cacheMove(uint64_t src, uint64_t dest, float mass) const
735
-{
736
-    unsigned ndx = 0;
737
-    #pragma omp critical(atomicMove)
738
-    {
739
-        ndx = mMoveCacheIndex++;
740
-    }
741
-    mMoveCache[ndx] = MoveOperation(src, dest, mass);
742
-}
743
-
744 141
 void AtomicDomain::resetCache(unsigned n)
745 142
 {
746 143
     mInsertCacheIndex = 0;
747 144
     mEraseCacheIndex = 0;
748
-    mMoveCacheIndex = 0;
749 145
     mInsertCache.resize(n);
750 146
     mEraseCache.resize(n);
751
-    mMoveCache.resize(n);
147
+}
148
+
149
+void AtomicDomain::erase(uint64_t pos)
150
+{
151
+    GAPS_ASSERT(size() > 0);
152
+    std::vector<Atom>::iterator it;
153
+    it = std::lower_bound(mAtoms.begin(), mAtoms.end(), pos, compareAtomLower);
154
+    mAtoms.erase(it);
155
+}
156
+
157
+void AtomicDomain::insert(uint64_t pos, float mass)
158
+{
159
+    std::vector<Atom>::iterator it;
160
+    it = std::lower_bound(mAtoms.begin(), mAtoms.end(), pos, compareAtomLower);
161
+    mAtoms.insert(it, Atom(pos, mass));
752 162
 }
753 163
 
754 164
 void AtomicDomain::flushCache()
755 165
 {
756 166
     for (unsigned i = 0; i < mEraseCacheIndex; ++i)
757 167
     {
758
-        mAtoms.erase(mEraseCache[i]);
168
+        erase(mEraseCache[i]);
759 169
     }
760 170
 
761 171
     for (unsigned i = 0; i < mInsertCacheIndex; ++i)
762 172
     {
763
-        mAtoms.insert(mInsertCache[i].pos, mInsertCache[i].mass);
764
-    }
765
-
766
-    for (unsigned i = 0; i < mMoveCacheIndex; ++i)
767
-    {
768
-        mAtoms.move(mMoveCache[i].src, mMoveCache[i].dest, mMoveCache[i].mass);
173
+        insert(mInsertCache[i].pos, mInsertCache[i].mass);
769 174
     }
770 175
 
771 176
     mInsertCache.clear();
772 177
     mEraseCache.clear();
773
-    mMoveCache.clear();
774 178
 }
775 179
 
776 180
 Archive& operator<<(Archive &ar, AtomicDomain &domain)
777 181
 {
778
-    ar << domain.mDomainLength << domain.mAtoms << domain.mRng;
182
+    ar << domain.mDomainLength << domain.mRng << domain.mAtoms.size();
183
+    
184
+    for (unsigned i = 0; i < domain.mAtoms.size(); ++i)
185
+    {
186
+        ar << domain.mAtoms[i];
187
+    }
188
+    return ar;
779 189
 }
780 190
 
781 191
 Archive& operator>>(Archive &ar, AtomicDomain &domain)
782 192
 {
783
-    ar >> domain.mDomainLength >> domain.mAtoms << domain.mRng;
193
+    Atom temp;
194
+    unsigned size = 0;
195
+    ar >> domain.mDomainLength >> domain.mRng >> size;
196
+    for (unsigned i = 0; i < size; ++i)
197
+    {
198
+        ar >> temp;
199
+        domain.insert(temp.pos, temp.mass);
200
+    }
201
+    return ar;
784 202
 }
... ...
@@ -5,7 +5,6 @@
5 5
 #include "math/Random.h"
6 6
 
7 7
 #include <vector>
8
-#include <set>
9 8
 
10 9
 struct Atom
11 10
 {
... ...
@@ -33,147 +32,12 @@ struct AtomNeighborhood
33 32
     bool hasRight();
34 33
 };
35 34
 
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
-    Atom* getLeft(unsigned index);
59
-    Atom* getRight(unsigned index);
60
-
61
-#ifndef GAPS_INTERNAL_TESTS
62
-private:
63
-#endif
64
-
65
-    Atom mBucket[2];
66
-    uint32_t mSize;
67
-    
68
-    AtomBucket *mOverflow;
69
-    AtomBucket *mPrev;
70
-    AtomBucket *mNext;
71
-
72
-    void eraseFront();
73
-    void connectAdjacent();
74
-
75
-    Atom* back();
76
-
77
-    AtomBucket& operator=(const AtomBucket& other); // prevent copies
78
-    
79
-    friend Archive& operator<<(Archive& ar, AtomBucket &b);
80
-    friend Archive& operator>>(Archive& ar, AtomBucket &b);
81
-};
82
-
83
-struct HashMapIndex
84
-{
85
-    unsigned bucket;
86
-    unsigned index;
87
-
88
-    HashMapIndex(unsigned b, unsigned i) : bucket(b), index(i) {}
89
-};
90
-
91
-// hash map of chained AtomBuckets
92
-// data structure that holds atoms
93
-// accessing happens much more than insert/erase so more time is spent
94
-// up front (sorting) to save time on access
95
-// note that atoms can never occupy position 0
96
-class AtomHashMap
97
-{
98
-public:
99
-
100
-    // required that length is a multiple of expectedNumAtoms
101
-    AtomHashMap(unsigned expectedNumAtoms);
102
-    void setTotalLength(uint64_t length);
103
-
104
-    // TODO while moving the atoms cannot change the atom order, it can
105
-    // move an atom from one bucket to another - this moved atom must be
106
-    // the front atom moving to the back of the prev bucket or the back
107
-    // atom moving to the front of the next atom - effectively creating
108
-    // a reorder event - it may even need allocation
109
-
110
-    Atom* front();
111
-    Atom* randomAtom();
112
-    AtomNeighborhood randomAtomWithNeighbors();
113
-    AtomNeighborhood randomAtomWithRightNeighbor();
114
-
115
-    bool contains(uint64_t pos) const;
116
-    unsigned size() const;
117
-
118
-    void insert(uint64_t pos, float mass);
119
-    void erase(uint64_t pos);
120
-    void move(uint64_t src, uint64_t dest, float mass); // cannot reorder
121
-
122
-    Atom* getLeft(uint64_t pos);
123
-    Atom* getRight(uint64_t pos);
124
-
125
-    bool hasLeft(uint64_t pos);
126
-    bool hasRight(uint64_t pos);
127
-
128
-    void updateMass(uint64_t pos, float newMass);
129
-
130
-#ifndef GAPS_INTERNAL_TESTS
131
-private:
132
-#endif
133
-
134
-    unsigned mSize;
135
-    unsigned mWhichLongestBucket;
136
-    unsigned mLongestBucketSize;
137
-
138
-    unsigned mExpectedNumAtoms;
139
-
140
-    uint64_t mBucketLength;
141
-    uint64_t mTotalLength;
142
-
143
-    std::vector<AtomBucket> mHashMap;
144
-    std::set<unsigned> mFullBuckets; // TODO use IntFixedHashSet
145
-
146
-    // random number generator
147
-    mutable GapsRng mRng;
148
-
149
-    unsigned hash(uint64_t pos) const;
150
-    HashMapIndex getRandomIndex() const;
151
-
152
-    friend Archive& operator<<(Archive& ar, AtomHashMap &h);
153
-    friend Archive& operator>>(Archive& ar, AtomHashMap &h);
154
-};
155
-
156
-struct MoveOperation
157
-{
158
-    uint64_t src;
159
-    uint64_t dest;
160
-    float mass;
161
-
162
-    MoveOperation() : src(0), dest(0), mass(0.f) {}
163
-
164
-    MoveOperation(uint64_t s, uint64_t d, float m) :
165
-        src(s), dest(d), mass(m)
166
-    {}
167
-};
168
-
169 35
 class AtomicDomain
170 36
 {
171 37
 public:
172 38
 
173 39
     AtomicDomain(unsigned nBins);
174 40
 
175
-    void setDomainSize(uint64_t size);
176
-
177 41
     // access atoms
178 42
     Atom* front();
179 43
     Atom* randomAtom();
... ...
@@ -186,7 +50,6 @@ public:
186 50
     // modify domain
187 51
     void cacheInsert(uint64_t pos, float mass) const;
188 52
     void cacheErase(uint64_t pos) const;
189
-    void cacheMove(uint64_t src, uint64_t dest, float mass) const;
190 53
     void flushCache();
191 54
     void resetCache(unsigned n);
192 55
 
... ...
@@ -194,27 +57,30 @@ public:
194 57
     friend Archive& operator<<(Archive &ar, AtomicDomain &domain);
195 58
     friend Archive& operator>>(Archive &ar, AtomicDomain &domain);
196 59
 
60
+#ifndef GAPS_INTERNAL_TESTS
197 61
 private:
62
+#endif
198 63
 
199 64
     // size of atomic domain to ensure all bins are equal length
200 65
     uint64_t mDomainLength;
201 66
 
202 67
     // domain storage, specialized hash map
203
-    mutable AtomHashMap mAtoms;
68
+    std::vector<Atom> mAtoms;
204 69
 
205 70
     // holds cache of operations
206 71
     mutable std::vector<Atom> mInsertCache;
207 72
     mutable std::vector<uint64_t> mEraseCache;
208
-    mutable std::vector<MoveOperation> mMoveCache;
209 73
 
210 74
     // current index in the operation cache
211 75
     mutable unsigned mInsertCacheIndex;
212 76
     mutable unsigned mEraseCacheIndex;
213
-    mutable unsigned mMoveCacheIndex;
214 77
 
215 78
     // random number generator
216 79
     mutable GapsRng mRng;
217 80
 
81
+    void erase(uint64_t pos);
82
+    void insert(uint64_t pos, float mass);
83
+
218 84
     // serialization
219 85
     friend Archive& operator<<(Archive &ar, AtomicDomain &domain);
220 86
     friend Archive& operator>>(Archive &ar, AtomicDomain &domain);
... ...
@@ -12,7 +12,6 @@
12 12
 #include <algorithm>
13 13
 
14 14
 // forward declarations needed for friend classes/functions
15
-
16 15
 class AmplitudeGibbsSampler;
17 16
 class PatternGibbsSampler;
18 17
 class GapsStatistics;
... ...
@@ -62,9 +61,6 @@ protected:
62 61
 
63 62
     void processProposal(AtomicProposal *prop);
64 63
 
65
-    void addMass(uint64_t pos, float mass, unsigned row, unsigned col);
66
-    void removeMass(uint64_t pos, float mass, unsigned row, unsigned col);
67
-
68 64
     void birth(AtomicProposal *prop);
69 65
     void death(AtomicProposal *prop);
70 66
     void move(AtomicProposal *prop);
... ...
@@ -184,9 +180,7 @@ bool transposeData, unsigned nPatterns, bool partitionRows,
184 180
 const std::vector<unsigned> &indices)
185 181
     :
186 182
 GibbsSampler(data, transposeData, nPatterns, true, partitionRows, indices)
187
-{
188
-    mQueue.setDimensionSize(mBinSize, mNumCols);
189
-}
183
+{}
190 184
 
191 185
 template <class DataType>
192 186
 PatternGibbsSampler::PatternGibbsSampler(const DataType &data,
... ...
@@ -194,9 +188,7 @@ bool transposeData, unsigned nPatterns, bool partitionRows,
194 188
 const std::vector<unsigned> &indices)
195 189
     :
196 190
 GibbsSampler(data, transposeData, nPatterns, false, partitionRows, indices)
197
-{
198
-    mQueue.setDimensionSize(mBinSize, mNumRows);
199
-}
191
+{}
200 192
 
201 193
 template <class T, class MatA, class MatB>
202 194
 template <class DataType>
... ...
@@ -208,7 +200,7 @@ mDMatrix(data, transposeData, partitionRows, indices),
208 200
 mSMatrix(mDMatrix.pmax(0.1f, 0.1f)),
209 201
 mAPMatrix(mDMatrix.nRow(), mDMatrix.nCol()),
210 202
 mMatrix(amp ? mDMatrix.nRow() : nPatterns, amp ? nPatterns : mDMatrix.nCol()),
211
-mOtherMatrix(NULL), mQueue(mMatrix.nRow() * mMatrix.nCol()),
203
+mOtherMatrix(NULL), mQueue(amp ? mDMatrix.nRow() : mDMatrix.nCol(), nPatterns),
212 204
 mDomain(mMatrix.nRow() * mMatrix.nCol()), mLambda(0.f),
213 205
 mMaxGibbsMass(100.f), mAnnealingTemp(1.f), mNumRows(mMatrix.nRow()),
214 206
 mNumCols(mMatrix.nCol()), mAvgQueue(0.f), mNumQueues(0.f)
... ...
@@ -216,8 +208,6 @@ mNumCols(mMatrix.nCol()), mAvgQueue(0.f), mNumQueues(0.f)
216 208
     // calculate atomic domain size
217 209
     mBinSize = std::numeric_limits<uint64_t>::max()
218 210
         / static_cast<uint64_t>(mNumRows * mNumCols);
219
-    mQueue.setDomainSize(mBinSize * mNumRows * mNumCols);
220
-    mDomain.setDomainSize(mBinSize * mNumRows * mNumCols);
221 211
 
222 212
     // default sparsity parameters
223 213
     setSparsity(0.01, false);
... ...
@@ -270,7 +260,8 @@ void GibbsSampler<T, MatA, MatB>::update(unsigned nSteps, unsigned nCores)
270 260
     while (n < nSteps)
271 261
     {
272 262
         // populate queue, prepare domain for this queue
273
-        mQueue.populate(mDomain, nSteps - n);
263
+        //mQueue.populate(mDomain, nSteps - n);
264
+        mQueue.populate(mDomain, 1);
274 265
         mDomain.resetCache(mQueue.size());
275 266
         n += mQueue.size();
276 267
         
... ...
@@ -287,6 +278,7 @@ void GibbsSampler<T, MatA, MatB>::update(unsigned nSteps, unsigned nCores)
287 278
         {
288 279
             processProposal(&mQueue[i]);
289 280
         }
281
+
290 282
         mDomain.flushCache();
291 283
         mQueue.clear();
292 284
     }
... ...
@@ -377,22 +369,6 @@ void GibbsSampler<T, MatA, MatB>::processProposal(AtomicProposal *prop)
377 369
     }
378 370
 }
379 371
 
380
-template <class T, class MatA, class MatB>
381
-void GibbsSampler<T, MatA, MatB>::addMass(uint64_t pos, float mass, unsigned row, unsigned col)
382
-{
383
-    mDomain.cacheInsert(pos, mass);
384
-    mMatrix(row, col) += mass;
385
-    impl()->updateAPMatrix(row, col, mass);
386
-}
387
-
388
-template <class T, class MatA, class MatB>
389
-void GibbsSampler<T, MatA, MatB>::removeMass(uint64_t pos, float mass, unsigned row, unsigned col)
390
-{
391
-    mDomain.cacheErase(pos);
392
-    mMatrix(row, col) = gaps::max(mMatrix(row, col) - mass, 0.f);
393
-    impl()->updateAPMatrix(row, col, -mass);
394
-}
395
-
396 372
 // add an atom at pos, calculate mass either with an exponential distribution
397 373
 // or with the gibbs mass calculation
398 374
 template <class T, class MatA, class MatB>
... ...
@@ -477,20 +453,20 @@ void GibbsSampler<T, MatA, MatB>::move(AtomicProposal *prop)
477 453
     unsigned r2 = impl()->getRow(prop->moveDest);
478 454
     unsigned c2 = impl()->getCol(prop->moveDest);
479 455
 
480
-    if (r1 != r2 || c1 != c2) // automatically reject if change in same bin
456
+    GAPS_ASSERT(r1 != r2 || c1 != c2); // same bin handled during proposal
457
+
458
+    float deltaLL = impl()->computeDeltaLL(r1, c1, -1.f * prop->atom1->mass,
459
+        r2, c2, prop->atom1->mass);
460
+    if (deltaLL * mAnnealingTemp > std::log(prop->rng.uniform()))
481 461
     {
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()))
485
-        {
486
-            mDomain.cacheMove(prop->atom1->pos, prop->moveDest, prop->atom1->mass);
462
+        //mDomain.cacheMove(prop->atom1->pos, prop->moveDest, prop->atom1->mass);
463
+        prop->atom1->pos = prop->moveDest; // temporary for vector domain
487 464
 
488
-            mMatrix(r1, c1) = gaps::max(mMatrix(r1, c1) - prop->atom1->mass, 0.f);
489
-            mMatrix(r2, c2) += prop->atom1->mass;
465
+        mMatrix(r1, c1) = gaps::max(mMatrix(r1, c1) - prop->atom1->mass, 0.f);
466
+        mMatrix(r2, c2) += prop->atom1->mass;
490 467
 
491
-            impl()->updateAPMatrix(r1, c1, -1.f * prop->atom1->mass);
492
-            impl()->updateAPMatrix(r2, c2, prop->atom1->mass);
493
-        }
468
+        impl()->updateAPMatrix(r1, c1, -1.f * prop->atom1->mass);
469
+        impl()->updateAPMatrix(r2, c2, prop->atom1->mass);
494 470
     }
495 471
 }
496 472
 
... ...
@@ -505,45 +481,46 @@ void GibbsSampler<T, MatA, MatB>::exchange(AtomicProposal *prop)
505 481
     unsigned r2 = impl()->getRow(prop->atom2->pos);
506 482
     unsigned c2 = impl()->getCol(prop->atom2->pos);
507 483
 
508
-    if (r1 != r2 || c1 != c2) // automatically reject if change in same bin
509
-    {
510
-        float m1 = prop->atom1->mass;
511
-        float m2 = prop->atom2->mass;
484
+    GAPS_ASSERT(r1 != r2 || c1 != c2); // same bin handled during proposal
512 485
 
513
-        if (impl()->canUseGibbs(r1, c1, r2, c2))
486
+    float m1 = prop->atom1->mass;
487
+    float m2 = prop->atom2->mass;
488
+
489
+    if (impl()->canUseGibbs(r1, c1, r2, c2))
490
+    {
491
+        AlphaParameters alpha = impl()->alphaParameters(r1, c1, r2, c2);
492
+        std::pair<float, bool> gMass = gibbsMass(alpha, m1, m2, &(prop->rng));
493
+        if (gMass.second)
514 494
         {
515
-            AlphaParameters alpha = impl()->alphaParameters(r1, c1, r2, c2);
516
-            std::pair<float, bool> gMass = gibbsMass(alpha, m1, m2, &(prop->rng));
517
-            if (gMass.second)
518
-            {
519
-                acceptExchange(prop, gMass.first, r1, c1, r2, c2);
520
-                return;
521
-            }
495
+            acceptExchange(prop, gMass.first, r1, c1, r2, c2);
496
+            return;
522 497
         }
498
+    }
523 499
 
524
-        float pUpper = gaps::p_gamma(m1 + m2, 2.f, 1.f / mLambda);
525
-        float newMass = prop->rng.inverseGammaSample(0.f, pUpper, 2.f, 1.f / mLambda);
500
+    float pUpper = gaps::p_gamma(m1 + m2, 2.f, 1.f / mLambda);
501
+    float newMass = prop->rng.inverseGammaSample(0.f, pUpper, 2.f, 1.f / mLambda);
526 502
 
527
-        float delta = m1 > m2 ? newMass - m1 : m2 - newMass; // change larger mass
528
-        float pOldMass = 2.f * newMass > m1 + m2 ? gaps::max(m1, m2) : gaps::min(m1, m2);
529
-    
530
-        float pNew = gaps::d_gamma(newMass, 2.f, 1.f / mLambda);
531
-        float pOld = gaps::d_gamma(pOldMass, 2.f, 1.f / mLambda);
503
+    float delta = m1 > m2 ? newMass - m1 : m2 - newMass; // change larger mass
504
+    float pOldMass = 2.f * newMass > m1 + m2 ? gaps::max(m1, m2) : gaps::min(m1, m2);
532 505
 
533
-        if (pOld == 0.f && pNew != 0.f) // special case
534
-        {
535
-            acceptExchange(prop, delta, r1, c1, r2, c2);
536
-            return;
537
-        }
538
-        float deltaLL = impl()->computeDeltaLL(r1, c1, delta, r2, c2, -delta);
539
-        float priorLL = (pOld == 0.f) ? 1.f : pOld / pNew;
540
-        float u = std::log(prop->rng.uniform() * priorLL);
541
-        if (u < deltaLL * mAnnealingTemp)
542
-        {
543
-            acceptExchange(prop, delta, r1, c1, r2, c2);
544
-            return;
545
-        }
506
+    float pNew = gaps::d_gamma(newMass, 2.f, 1.f / mLambda);
507
+    float pOld = gaps::d_gamma(pOldMass, 2.f, 1.f / mLambda);
508
+
509
+    if (pOld == 0.f && pNew != 0.f) // special case
510
+    {
511
+        acceptExchange(prop, delta, r1, c1, r2, c2);
512
+        return;
546 513
     }
514
+
515
+    float deltaLL = impl()->computeDeltaLL(r1, c1, delta, r2, c2, -delta);
516
+    float priorLL = (pOld == 0.f) ? 1.f : pOld / pNew;
517
+    float u = std::log(prop->rng.uniform() * priorLL);
518
+    if (u < deltaLL * mAnnealingTemp)
519
+    {
520
+        acceptExchange(prop, delta, r1, c1, r2, c2);
521
+        return;
522
+    }
523
+
547 524
     mQueue.rejectDeath();
548 525
 }
549 526
 
... ...
@@ -568,6 +545,12 @@ template <class T, class MatA, class MatB>
568 545
 void GibbsSampler<T, MatA, MatB>::acceptExchange(AtomicProposal *prop,
569 546
 float d1, unsigned r1, unsigned c1, unsigned r2, unsigned c2)
570 547
 {
548
+    // can make this assertion more formal if we customize p/q norm functions
549
+    // to handle tail situations a certain way
550
+    GAPS_ASSERT(prop->atom1->mass + d1 > gaps::epsilon);
551
+    GAPS_ASSERT(prop->atom2->mass - d1 > gaps::epsilon);
552
+    //d1 = gaps::max(-1.f * (prop->atom1->mass - epsilon), d1);
553
+
571 554
     float d2 = -1.f * d1;
572 555
     bool b1 = updateAtomMass(prop->atom1, d1);
573 556
     bool b2 = updateAtomMass(prop->atom2, d2);
... ...
@@ -2,20 +2,14 @@
2 2
 #include "ProposalQueue.h"
3 3
 #include "math/Random.h"
4 4
 
5
-void ProposalQueue::setNumBins(unsigned nBins)
5
+ProposalQueue::ProposalQueue(unsigned primaryDimSize, unsigned secondaryDimSize)
6
+    :
7
+mMinAtoms(0), mMaxAtoms(0), mNumBins(primaryDimSize * secondaryDimSize),
8
+mBinLength(std::numeric_limits<uint64_t>::max() / mNumBins),
9
+mSecondaryDimLength(mBinLength * secondaryDimSize),
10
+mDomainLength(mBinLength * mNumBins), mSecondaryDimSize(secondaryDimSize)
6 11
 {
7
-    mNumBins = nBins;
8
-}
9
-
10
-void ProposalQueue::setDomainSize(uint64_t size)
11
-{
12
-    mDomainSize = size;
13
-}
14
-
15
-void ProposalQueue::setDimensionSize(uint64_t binSize, uint64_t dimLength)
16
-{
17
-    mDimensionSize = binSize * dimLength;
18
-    mUsedIndices.setDimensionSize(mNumBins / dimLength);
12
+    mUsedIndices.setDimensionSize(primaryDimSize);
19 13
 }
20 14
 
21 15
 void ProposalQueue::setAlpha(float alpha)
... ...
@@ -81,7 +75,7 @@ void ProposalQueue::rejectBirth()
81 75
 
82 76
 float ProposalQueue::deathProb(uint64_t nAtoms) const
83 77
 {
84
-    double size = static_cast<double>(mDomainSize);
78
+    double size = static_cast<double>(mDomainLength);
85 79
     double term1 = (size - static_cast<double>(nAtoms)) / size;
86 80
     double term2 = mAlpha * static_cast<double>(mNumBins) * term1;
87 81
     return static_cast<double>(nAtoms) / (static_cast<double>(nAtoms) + term2);
... ...
@@ -120,21 +114,33 @@ bool ProposalQueue::makeProposal(AtomicDomain &domain)
120 114
         {
121 115
             return death(domain);
122 116
         }
123
-        return false;
117
+        return false; // can't determine B/D since range is too wide
124 118
     }
125 119
     return (mU1 < 0.75f || mMaxAtoms < 2) ? move(domain) : exchange(domain);
126 120
 }
127 121
     
122
+unsigned ProposalQueue::primaryIndex(uint64_t pos) const
123
+{
124
+    return pos / mSecondaryDimLength;
125
+}
126
+
127
+unsigned ProposalQueue::secondaryIndex(uint64_t pos) const
128
+{
129
+    return (pos / mBinLength) % mSecondaryDimSize;
130
+}
131
+
132
+// TODO add atoms with empty mass? fill in mass in gibbssampler?
133
+// inserting invalidates previous pointers
128 134
 bool ProposalQueue::birth(AtomicDomain &domain)
129 135
 {
130 136
     uint64_t pos = domain.randomFreePosition();
131
-    if (mUsedIndices.count(pos / mDimensionSize))
137
+    if (mUsedIndices.contains(primaryIndex(pos)))
132 138
     {
133 139
         return false; // matrix conflict - can't compute gibbs mass
134 140
     }
135 141
 
136 142
     mQueue.push_back(AtomicProposal('B', pos));
137
-    mUsedIndices.insert(pos / mDimensionSize);
143
+    mUsedIndices.insert(primaryIndex(pos));
138 144
     mUsedPositions.insert(pos);
139 145
     ++mMaxAtoms;
140 146
     return true;
... ...
@@ -143,13 +149,13 @@ bool ProposalQueue::birth(AtomicDomain &domain)
143 149
 bool ProposalQueue::death(AtomicDomain &domain)
144 150
 {
145 151
     Atom* a = domain.randomAtom();
146
-    if (mUsedIndices.count(a->pos / mDimensionSize))
152
+    if (mUsedIndices.contains(primaryIndex(a->pos)))
147 153
     {
148 154
         return false; // matrix conflict - can't compute gibbs mass or deltaLL
149 155
     }
150 156
 
151 157
     mQueue.push_back(AtomicProposal('D', a));
152
-    mUsedIndices.insert(a->pos / mDimensionSize);
158
+    mUsedIndices.insert(primaryIndex(a->pos));
153 159
     mUsedPositions.insert(a->pos);
154 160
     --mMinAtoms;
155 161
     return true;
... ...
@@ -158,23 +164,32 @@ bool ProposalQueue::death(AtomicDomain &domain)
158 164
 bool ProposalQueue::move(AtomicDomain &domain)
159 165
 {
160 166
     AtomNeighborhood hood = domain.randomAtomWithNeighbors();
161
-    uint64_t lbound = hood.hasLeft() ? hood.left->pos : 1;
162
-    uint64_t rbound = hood.hasRight() ? hood.right->pos : mDomainSize;
167
+    uint64_t lbound = hood.hasLeft() ? hood.left->pos : 0;
168
+    uint64_t rbound = hood.hasRight() ? hood.right->pos : mDomainLength;
163 169
 
164 170
     if (!mUsedPositions.isEmptyInterval(lbound, rbound))
165 171
     {
166
-        return false;
172
+        return false; // atomic conflict - don't know neighbors
173
+    }
174
+
175
+    uint64_t newLocation = mRng.uniform64(lbound + 1, rbound - 1);
176
+
177
+    if (primaryIndex(hood.center->pos) == primaryIndex(newLocation)
178
+    && secondaryIndex(hood.center->pos) == secondaryIndex(newLocation))
179
+    {
180
+        hood.center->pos = newLocation; // automatically accept moves in same bin
181
+        return true;
167 182
     }
168 183
 
169
-    uint64_t newLocation = mRng.uniform64(lbound, rbound - 1);
170
-    if (mUsedIndices.count(hood.center->pos / mDimensionSize) || mUsedIndices.count(newLocation / mDimensionSize))
184
+    if (mUsedIndices.contains(primaryIndex(hood.center->pos))
185
+    || mUsedIndices.contains(primaryIndex(newLocation)))
171 186
     {
172 187
         return false; // matrix conflict - can't compute deltaLL
173 188
     }
174 189
 
175 190
     mQueue.push_back(AtomicProposal('M', hood.center, newLocation));
176
-    mUsedIndices.insert(hood.center->pos / mDimensionSize);
177
-    mUsedIndices.insert(newLocation / mDimensionSize);
191
+    mUsedIndices.insert(primaryIndex(hood.center->pos));
192
+    mUsedIndices.insert(primaryIndex(newLocation));
178 193
     mUsedPositions.insert(hood.center->pos);
179 194
     mUsedPositions.insert(newLocation);
180 195
     return true;
... ...
@@ -190,30 +205,37 @@ bool ProposalQueue::exchange(AtomicDomain &domain)
190 205
     {
191 206
         if (!mUsedPositions.isEmptyInterval(a1->pos, a2->pos))
192 207
         {
193
-            return false;
208
+            return false; // atomic conflict - don't know right neighbor
194 209
         }
195 210
     }
196 211
     else // exchange with first atom
197 212
     {
198
-        if (!mUsedPositions.isEmptyInterval(a1->pos, mDomainSize))
213
+        if (!mUsedPositions.isEmptyInterval(a1->pos, mDomainLength))
199 214
         {
200
-            return false;
215
+            return false; // atomic conflict - don't know right neighbor
201 216
         }
202 217
         
203 218
         if (!mUsedPositions.isEmptyInterval(0, domain.front()->pos))
204 219
         {
205
-            return false;
220
+            return false; // atomic conflict - don't know right neighbor
206 221
         }
207 222
     }
208 223
 
209
-    if (mUsedIndices.count(a1->pos / mDimensionSize) || mUsedIndices.count(a2->pos / mDimensionSize))
224
+    if (primaryIndex(a1->pos) == primaryIndex(a2->pos)
225
+    && secondaryIndex(a1->pos) == secondaryIndex(a2->pos))
226
+    {
227
+        return true; // TODO automatically accept exchanges in same bin
228
+    }
229
+
230
+    if (mUsedIndices.contains(primaryIndex(a1->pos))
231
+    || mUsedIndices.contains(primaryIndex(a2->pos)))
210 232
     {
211 233
         return false; // matrix conflict - can't compute gibbs mass or deltaLL
212 234
     }
213 235
 
214 236
     mQueue.push_back(AtomicProposal('E', a1, a2));
215
-    mUsedIndices.insert(a1->pos / mDimensionSize);
216
-    mUsedIndices.insert(a2->pos / mDimensionSize);
237
+    mUsedIndices.insert(primaryIndex(a1->pos));
238
+    mUsedIndices.insert(primaryIndex(a2->pos));
217 239
     mUsedPositions.insert(a1->pos);
218 240
     mUsedPositions.insert(a2->pos);
219 241
     --mMinAtoms;
... ...
@@ -223,15 +245,17 @@ bool ProposalQueue::exchange(AtomicDomain &domain)
223 245
 Archive& operator<<(Archive &ar, ProposalQueue &queue)
224 246
 {
225 247
     ar << queue.mMinAtoms << queue.mMaxAtoms << queue.mNumBins
226
-        << queue.mDimensionSize << queue.mDomainSize << queue.mAlpha
227
-        << queue.mRng;
248
+        << queue.mBinLength << queue.mSecondaryDimLength << queue.mDomainLength
249
+        << queue.mSecondaryDimSize << queue.mAlpha << queue.mRng
250
+        << queue.mU1 << queue.mU2 << queue.mUseCachedRng;
228 251
     return ar;
229 252
 }
230 253
 
231 254
 Archive& operator>>(Archive &ar, ProposalQueue &queue)
232 255
 {
233 256
     ar >> queue.mMinAtoms >> queue.mMaxAtoms >> queue.mNumBins
234
-        >> queue.mDimensionSize >> queue.mDomainSize >> queue.mAlpha
235
-        >> queue.mRng;
257
+        >> queue.mBinLength >> queue.mSecondaryDimLength >> queue.mDomainLength
258
+        >> queue.mSecondaryDimSize >> queue.mAlpha >> queue.mRng
259
+        >> queue.mU1 >> queue.mU2 >> queue.mUseCachedRng;
236 260
     return ar;
237 261
 }
... ...
@@ -6,10 +6,9 @@
6 6
 #include "data_structures/EfficientSets.h"
7 7
 #include "math/Random.h"
8 8
 
9
-#include <boost/unordered_set.hpp>
10
-
11 9
 #include <cstddef>
12 10
 #include <stdint.h>
11
+#include <vector>
13 12
 
14 13
 struct AtomicProposal
15 14
 {
... ...
@@ -45,6 +44,23 @@ struct AtomicProposal
45 44
 
46 45
 class ProposalQueue
47 46
 {
47
+public:
48
+
49
+    ProposalQueue(unsigned primaryDimSize, unsigned secondaryDimSize);
50
+    void setAlpha(float alpha);
51
+
52
+    // modify/access queue
53
+    void populate(AtomicDomain &domain, unsigned limit);
54
+    void clear();
55
+    unsigned size() const;
56
+    AtomicProposal& operator[](int n);
57
+
58
+    // update min/max atoms
59
+    void acceptDeath();
60
+    void rejectDeath();
61
+    void acceptBirth();
62
+    void rejectBirth();
63
+
48 64
 private:
49 65
 
50 66
     std::vector<AtomicProposal> mQueue; // not really a queue for now
... ...
@@ -55,17 +71,22 @@ private:
55 71
     uint64_t mMinAtoms;
56 72
     uint64_t mMaxAtoms;
57 73
 
58
-    uint64_t mNumBins;
59
-    uint64_t mDimensionSize; // rows of A, cols of P
60
-    uint64_t mDomainSize;
74
+    unsigned mNumBins; // number of matrix elements
75
+    uint64_t mBinLength; // atomic length of one bin
76
+    uint64_t mSecondaryDimLength; // atomic length of one row (col) for A (P)
77
+    uint64_t mDomainLength; // length of entire atomic domain
78
+    unsigned mSecondaryDimSize; // number of cols (rows) for A (P)
61 79
 
62 80
     float mAlpha;
63 81
 
64
-    bool mUseCachedRng;
82
+    mutable GapsRng mRng;
83
+
65 84
     float mU1;
66 85
     float mU2;
86
+    bool mUseCachedRng;
67 87
 
68
-    mutable GapsRng mRng;
88
+    unsigned primaryIndex(uint64_t pos) const;
89
+    unsigned secondaryIndex(uint64_t pos) const;
69 90
 
70 91
     float deathProb(uint64_t nAtoms) const;
71 92
     bool birth(AtomicDomain &domain);
... ...
@@ -75,31 +96,6 @@ private:
75 96
 
76 97
     bool makeProposal(AtomicDomain &domain);
77 98
 
78
-public:
79
-
80
-    ProposalQueue(unsigned nBins)
81
-        : mMinAtoms(0), mMaxAtoms(0), mNumBins(nBins), mDimensionSize(0),
82
-        mDomainSize(0), mAlpha(0.f), mUseCachedRng(false), mU1(0.f), mU2(0.f)
83
-    {}
84
-
85
-    // set parameters
86
-    void setNumBins(unsigned nBins);
87
-    void setDomainSize(uint64_t size);
88
-    void setAlpha(float alpha);
89
-    void setDimensionSize(uint64_t binSize, uint64_t dimLength);
90
-
91
-    // modify/access queue
92
-    void populate(AtomicDomain &domain, unsigned limit);
93
-    void clear();
94
-    unsigned size() const;
95
-    AtomicProposal& operator[](int n);
96
-
97
-    // update min/max atoms
98
-    void acceptDeath();
99
-    void rejectDeath();
100
-    void acceptBirth();
101
-    void rejectBirth();
102
-
103 99
     // serialization
104 100
     friend Archive& operator<<(Archive &ar, ProposalQueue &queue);
105 101
     friend Archive& operator>>(Archive &ar, ProposalQueue &queue);
... ...
@@ -2,6 +2,110 @@
2 2
 #include "../AtomicDomain.h"
3 3
 #include "../GapsPrint.h"
4 4
 
5
+TEST_CASE("AtomicDomain")
6
+{
7
+    SECTION("Construction")
8
+    {
9
+        AtomicDomain domain(10);
10
+    
11
+        REQUIRE(domain.size() == 0);
12
+    }
13
+
14
+    #ifdef GAPS_INTERNAL_TESTS
15
+    SECTION("Insert")
16
+    {
17
+        AtomicDomain domain(10);
18
+        
19
+        for (unsigned i = 0; i < 1000; ++i)
20
+        {
21
+            REQUIRE_NOTHROW(domain.insert(i, static_cast<float>(i)));
22
+            REQUIRE(domain.size() == i + 1);
23
+        }
24
+    }
25
+
26
+    SECTION("Erase")
27
+    {
28
+        AtomicDomain domain(10);
29
+        
30
+        for (unsigned i = 0; i < 1000; ++i)
31
+        {
32
+            domain.insert(i, static_cast<float>(i));
33
+        }
34
+
35
+        unsigned counter = domain.size();
36
+        for (unsigned j = 0; j < 1000; j += 10)
37
+        {
38
+            REQUIRE_NOTHROW(domain.erase(j));
39
+            REQUIRE(domain.size() == --counter);
40
+        }
41
+    }
42
+    
43
+    SECTION("Random Free Position")
44
+    {
45
+        AtomicDomain domain(10);
46
+
47
+        for (unsigned i = 0; i < 1000; ++i)
48
+        {
49
+            domain.insert(i, static_cast<float>(i));
50
+            REQUIRE_NOTHROW(domain.randomFreePosition());
51
+        }
52
+    }
53
+
54
+    SECTION("Random Atom")
55
+    {
56
+        AtomicDomain domain(10);
57
+
58
+        for (unsigned i = 0; i < 1000; ++i)
59
+        {
60
+            domain.insert(i, static_cast<float>(i));
61
+            
62
+            // single random atom
63
+            REQUIRE(domain.randomAtom()->pos < i + 1);
64
+
65
+            // random atom for exchange
66
+            AtomNeighborhood hood = domain.randomAtomWithRightNeighbor();
67
+            REQUIRE(hood.center->pos < i + 1);
68
+
69
+            REQUIRE(!hood.hasLeft());
70
+            if (hood.center->pos == i)
71
+            {
72
+                REQUIRE(!hood.hasRight());
73
+            }
74
+            else
75
+            {
76
+                REQUIRE(hood.hasRight());
77
+                REQUIRE(hood.right->pos == hood.center->pos + 1);
78
+            }
79
+
80
+            // random atom for move
81
+            hood = domain.randomAtomWithNeighbors();
82
+            REQUIRE(hood.center->pos < i + 1);
83
+
84
+            if (hood.center->pos == 0)
85
+            {
86
+                REQUIRE(!hood.hasLeft());
87
+            }
88
+            else
89
+            {
90
+                REQUIRE(hood.left->pos == hood.center->pos - 1);
91
+            }
92
+
93
+            if (hood.center->pos == i)
94
+            {
95
+                REQUIRE(!hood.hasRight());
96
+            }
97
+            else
98
+            {
99
+                REQUIRE(hood.hasRight());
100
+                REQUIRE(hood.right->pos == hood.center->pos + 1);
101
+            }
102
+        }
103
+    }
104
+    #endif
105
+}
106
+
107
+#if 0
108
+
5 109
 // used to create aligned buckets for testing
6 110
 AtomBucket* getAlignedBucket()
7 111
 {
... ...
@@ -652,4 +756,6 @@ TEST_CASE("AtomHashMap")
652 756
         REQUIRE(sum > 0.95f * 57000.f);
653 757
         REQUIRE(sum < 1.05f * 57000.f);
654 758
     }
655
-}
656 759
\ No newline at end of file
760
+}
761
+
762
+#endif
657 763
\ No newline at end of file
... ...
@@ -6,16 +6,16 @@ TEST_CASE("Test IntFixedHashSet")
6 6
     IntFixedHashSet set;
7 7
 
8 8
     set.setDimensionSize(1000);
9
-    REQUIRE(!set.count(123));
9
+    REQUIRE(!set.contains(123));
10 10
     
11 11
     set.insert(123);
12
-    REQUIRE(set.count(123));
12
+    REQUIRE(set.contains(123));
13 13
 
14 14
     set.clear();
15
-    REQUIRE(!set.count(123));
15
+    REQUIRE(!set.contains(123));
16 16
 
17 17
     set.insert(123);
18
-    REQUIRE(set.count(123));
18
+    REQUIRE(set.contains(123));
19 19
 }
20 20
 
21 21
 TEST_CASE("Test IntDenseOrderedSet")
... ...
@@ -16,10 +16,10 @@ public:
16 16
 
17 17
     IntFixedHashSet() : mCurrentKey(1) {}
18 18
 
19
-    void setDimensionSize(uint64_t size) {mSet.resize(size, 0);}
19
+    void setDimensionSize(unsigned size) {mSet.resize(size, 0);}
20 20
     void clear() {++mCurrentKey;}
21
-    bool count(uint64_t n) {return mSet[n] == mCurrentKey;}
22
-    void insert(uint64_t n) {mSet[n] = mCurrentKey;}
21
+    bool contains(unsigned n) {return mSet[n] == mCurrentKey;}
22
+    void insert(unsigned n) {mSet[n] = mCurrentKey;}
23 23
 };
24 24
 
25 25
 // TODO have sorted vector with at least some % of holes
... ...
@@ -10,13 +10,6 @@
10 10
 #include <boost/math/distributions/gamma.hpp>
11 11
 #include <boost/math/distributions/normal.hpp>
12 12
 
13
-#include <boost/random/exponential_distribution.hpp>
14
-#include <boost/random/mersenne_twister.hpp>
15
-#include <boost/random/normal_distribution.hpp>
16
-#include <boost/random/poisson_distribution.hpp>
17
-#include <boost/random/uniform_01.hpp>
18
-#include <boost/random/uniform_real_distribution.hpp>
19
-
20 13
 #include <algorithm>
21 14
 #include <stdint.h>
22 15