Browse code

debugging

sherman5 authored on 17/04/2018 18:45:37
Showing15 changed files

... ...
@@ -95,210 +95,125 @@ bool gaps::algo::isColZero(const ColMatrix &mat, unsigned col)
95 95
     return gaps::algo::sum(mat.getCol(col)) == 0;
96 96
 }
97 97
 
98
-/*
99
-// horribly slow, don't call often
100
-void gaps::algo::matrixMultiplication(TwoWayMatrix &C, const ColMatrix &A,
101
-const RowMatrix &B)
102
-{
103
-    for (unsigned i = 0; i < C.nRow(); ++i)
104
-    {
105
-        for (unsigned j = 0; j < C.nCol(); ++j)
106
-        {
107
-            float sum = 0.0;
108
-            for (unsigned k = 0; k < A.nCol(); ++k)
109
-            {
110
-                sum += A(i,k) * B(k,j);
111
-            }
112
-            C.set(i, j, sum);
113
-        }
114
-    }
115
-}
116
-
117
-static float deltaLL_comp(unsigned size, const float *D, const float *S,
118
-const float *AP, const float *other, float delta)
98
+AlphaParameters gaps::algo::alphaParameters(unsigned size, const float *D,
99
+const float *S, const float *AP, const float *mat)
119 100
 {
120
-    const gaps::simd::packedFloat pDelta = delta, two = 2.f;
121
-    gaps::simd::packedFloat d, pOther, pD, pAP, pS, partialSum = 0.f;
101
+    gaps::simd::packedFloat ratio, pMat, pD, pAP, pS;
102
+    gaps::simd::packedFloat partialS = 0.f, partialSU = 0.f;
122 103
     gaps::simd::Index i = 0;
123 104
     for (; i <= size - i.increment(); ++i)
124 105
     {   
125
-        pOther.load(other + i);
106
+        pMat.load(mat + i);
126 107
         pD.load(D + i);
127 108
         pAP.load(AP + i);
128 109
         pS.load(S + i);
129
-        d = pDelta * pOther;
130
-        partialSum += (d * (two * (pD - pAP) - d)) / (two * pS * pS);
110
+        ratio = pMat / pS;
111
+        partialS += ratio * ratio;
112
+        partialSU += ratio * (pD - pAP) / pS;
131 113
     }
132
-    float fd, delLL = partialSum.scalar();
114
+    float fratio, s = partialS.scalar(), su = partialSU.scalar();
133 115
     for (unsigned j = i.value(); j < size; ++j)
134 116
     {
135
-        fd = delta * other[j];
136
-        delLL += (fd * (2.f * (D[j] - AP[j]) - fd)) / (2.f * S[j] * S[j]);
117
+        fratio = mat[j] / S[j];
118
+        s += fratio * fratio;
119
+        su += fratio * (D[j] - AP[j]) / S[j];
137 120
     }
138
-    return delLL;
121
+    return AlphaParameters(s,su);
139 122
 }
140 123
 
141
-static float deltaLL_comp(unsigned size, const float *D, const float *S,
142
-const float *AP, const float *other1, float delta1, const float *other2,
143
-float delta2)
124
+AlphaParameters gaps::algo::alphaParameters(unsigned size, const float *D,
125
+const float *S, const float *AP, const float *mat1, const float *mat2)
144 126
 {
145
-    const gaps::simd::packedFloat pDelta1 = delta1, pDelta2 = delta2, two = 2.f;
146
-    gaps::simd::packedFloat d, pOther1, pOther2, pD, pAP, pS, partialSum = 0.f;
127
+    gaps::simd::packedFloat ratio, pMat1, pMat2, pD, pAP, pS;
128
+    gaps::simd::packedFloat partialS = 0.f, partialSU = 0.f;
147 129
     gaps::simd::Index i = 0;
148 130
     for (; i <= size - i.increment(); ++i)
149 131
     {   
150
-        pOther1.load(other1 + i);
151
-        pOther2.load(other2 + i);
132
+        pMat1.load(mat1 + i);
133
+        pMat2.load(mat2 + i);
152 134
         pD.load(D + i);
153 135
         pAP.load(AP + i);
154 136
         pS.load(S + i);
155
-        d = pDelta1 * pOther1 + pDelta2 * pOther2;
156
-        partialSum += (d * (two * (pD - pAP) - d)) / (two * pS * pS);
137
+        ratio = (pMat1 - pMat2) / pS;
138
+        partialS += ratio * ratio;
139
+        partialSU += ratio * (pD - pAP) / pS;
157 140
     }
158
-    float fd, delLL = partialSum.scalar();
141
+    float fratio, s = partialS.scalar(), su = partialSU.scalar();
159 142
     for (unsigned j = i.value(); j < size; ++j)
160 143
     {
161
-        fd = delta1 * other1[j] + delta2 * other2[j];
162
-        delLL += (fd * (2.f * (D[j] - AP[j]) - fd)) / (2.f * S[j] * S[j]);
163
-    }
164
-    return delLL;
165
-}
166
-
167
-float gaps::algo::deltaLL(const MatrixChange &ch, const TwoWayMatrix &D,
168
-const TwoWayMatrix &S, const ColMatrix &A, const RowMatrix &P,
169
-const TwoWayMatrix &AP)
170
-{
171
-    if (ch.label == 'A' && ch.nChanges == 2 && ch.row1 == ch.row2)
172
-    {
173
-        return deltaLL_comp(D.nCol(), D.rowPtr(ch.row1), S.rowPtr(ch.row1),
174
-            AP.rowPtr(ch.row1), P.rowPtr(ch.col1), ch.delta1,
175
-            P.rowPtr(ch.col2), ch.delta2);
176
-    }
177
-    else if (ch.label == 'A') // single change, or two independent changes
178
-    {
179
-        float d1 = 0.f, d2 = 0.f;
180
-        d1 = deltaLL_comp(D.nCol(), D.rowPtr(ch.row1), S.rowPtr(ch.row1),
181
-            AP.rowPtr(ch.row1), P.rowPtr(ch.col1), ch.delta1);
182
-        if (ch.nChanges == 2)
183
-        {
184
-            d2 = deltaLL_comp(D.nCol(), D.rowPtr(ch.row2), S.rowPtr(ch.row2),
185
-                AP.rowPtr(ch.row2), P.rowPtr(ch.col2), ch.delta2);
186
-        }
187
-        return d1 + d2;
188
-    }
189
-    else if (ch.label == 'P' && ch.nChanges == 2 && ch.col1 == ch.col2)
190
-    {
191
-        return deltaLL_comp(D.nRow(), D.colPtr(ch.col1), S.colPtr(ch.col1),
192
-            AP.colPtr(ch.col1), A.colPtr(ch.row1), ch.delta1,
193
-            A.colPtr(ch.row2), ch.delta2);
194
-    }
195
-    else // single change, or two independent changes
196
-    {
197
-        float d1 = 0.f, d2 = 0.f;
198
-        d1 = deltaLL_comp(D.nRow(), D.colPtr(ch.col1), S.colPtr(ch.col1),
199
-            AP.colPtr(ch.col1), A.colPtr(ch.row1), ch.delta1);
200
-        if (ch.nChanges == 2)
201
-        {
202
-            d2 = deltaLL_comp(D.nRow(), D.colPtr(ch.col2), S.colPtr(ch.col2),
203
-                AP.colPtr(ch.col2), A.colPtr(ch.row2), ch.delta2);
204
-        }
205
-        return d1 + d2;
144
+        fratio = (mat1[j] - mat2[j]) / S[j];
145
+        s += fratio * fratio;
146
+        su += fratio * (D[j] - AP[j]) / S[j];
206 147
     }
148
+    return AlphaParameters(s,su);
207 149
 }
208 150
 
209
-// single change
210
-static AlphaParameters alphaParameters_comp(unsigned size, const float *D,
211
-const float *S, const float *AP, const float *other)
151
+float gaps::algo::deltaLL(unsigned size, const float *D, const float *S,
152
+const float *AP, const float *mat, float delta)
212 153
 {
213
-    gaps::simd::packedFloat ratio, pOther, pD, pAP, pS;
214
-    gaps::simd::packedFloat partialS = 0.f, partialSU = 0.f;
154
+    const gaps::simd::packedFloat pDelta = delta, two = 2.f;
155
+    gaps::simd::packedFloat d, pMat, pD, pAP, pS, partialSum = 0.f;
215 156
     gaps::simd::Index i = 0;
216 157
     for (; i <= size - i.increment(); ++i)
217 158
     {   
218
-        pOther.load(other + i);
159
+        pMat.load(mat + i);
219 160
         pD.load(D + i);
220 161
         pAP.load(AP + i);
221 162
         pS.load(S + i);
222
-        ratio = pOther / pS;
223
-        partialS += ratio * ratio;
224
-        partialSU += ratio * (pD - pAP) / pS;
163
+        d = pDelta * pMat;
164
+        partialSum += (d * (two * (pD - pAP) - d)) / (two * pS * pS);
225 165
     }
226
-    float fratio, s = partialS.scalar(), su = partialSU.scalar();
166
+    float fd, delLL = partialSum.scalar();
227 167
     for (unsigned j = i.value(); j < size; ++j)
228 168
     {
229
-        fratio = other[j] / S[j];
230
-        s += fratio * fratio;
231
-        su += fratio * (D[j] - AP[j]) / S[j];
169
+        fd = delta * mat[j];
170
+        delLL += (fd * (2.f * (D[j] - AP[j]) - fd)) / (2.f * S[j] * S[j]);
232 171
     }
233
-    return AlphaParameters(s,su);
172
+    return delLL;
234 173
 }
235 174
 
236
-// two dependent changes
237
-static AlphaParameters alphaParameters_comp(unsigned size, const float *D,
238
-const float *S, const float *AP, const float *other1, const float *other2)
175
+float gaps::algo::deltaLL(unsigned size, const float *D, const float *S,
176
+const float *AP, const float *mat1, float delta1, const float *mat2,
177
+float delta2)
239 178
 {
240
-    gaps::simd::packedFloat ratio, pOther1, pOther2, pD, pAP, pS;
241
-    gaps::simd::packedFloat partialS = 0.f, partialSU = 0.f;
179
+    const gaps::simd::packedFloat pDelta1 = delta1, pDelta2 = delta2, two = 2.f;
180
+    gaps::simd::packedFloat d, pMat1, pMat2, pD, pAP, pS, partialSum = 0.f;
242 181
     gaps::simd::Index i = 0;
243 182
     for (; i <= size - i.increment(); ++i)
244 183
     {   
245
-        pOther1.load(other1 + i);
246
-        pOther2.load(other2 + i);
184
+        pMat1.load(mat1 + i);
185
+        pMat2.load(mat2 + i);
247 186
         pD.load(D + i);
248 187
         pAP.load(AP + i);
249 188
         pS.load(S + i);
250
-        ratio = (pOther1 - pOther2) / pS;
251
-        partialS += ratio * ratio;
252
-        partialSU += ratio * (pD - pAP) / pS;
189
+        d = pDelta1 * pMat1 + pDelta2 * pMat2;
190
+        partialSum += (d * (two * (pD - pAP) - d)) / (two * pS * pS);
253 191
     }
254
-    float fratio, s = partialS.scalar(), su = partialSU.scalar();
192
+    float fd, delLL = partialSum.scalar();
255 193
     for (unsigned j = i.value(); j < size; ++j)
256 194
     {
257
-        fratio = (other1[j] - other2[j]) / S[j];
258
-        s += fratio * fratio;
259
-        su += fratio * (D[j] - AP[j]) / S[j];
195
+        fd = delta1 * mat1[j] + delta2 * mat2[j];
196
+        delLL += (fd * (2.f * (D[j] - AP[j]) - fd)) / (2.f * S[j] * S[j]);
260 197
     }
261
-    return AlphaParameters(s,su);
198
+    return delLL;
262 199
 }
263 200
 
264
-// should these always be positive? could explain ordering of atoms in exchange
265
-AlphaParameters gaps::algo::alphaParameters(const MatrixChange &ch,
266
-const TwoWayMatrix &D, const TwoWayMatrix &S, const ColMatrix &A,
267
-const RowMatrix &P, const TwoWayMatrix &AP)
201
+/*
202
+// horribly slow, don't call often
203
+void gaps::algo::matrixMultiplication(TwoWayMatrix &C, const ColMatrix &A,
204
+const RowMatrix &B)
268 205
 {
269
-    if (ch.label == 'A' && ch.nChanges == 2 && ch.row1 == ch.row2)
270
-    {
271
-        return alphaParameters_comp(D.nCol(), D.rowPtr(ch.row1), S.rowPtr(ch.row1),
272
-            AP.rowPtr(ch.row1), P.rowPtr(ch.col1), P.rowPtr(ch.col2));
273
-    }
274
-    else if (ch.label == 'A') // single change, or two independent changes
275
-    {
276
-        AlphaParameters a1(0.f, 0.f), a2(0.f, 0.f);
277
-        a1 = alphaParameters_comp(D.nCol(), D.rowPtr(ch.row1), S.rowPtr(ch.row1),
278
-            AP.rowPtr(ch.row1), P.rowPtr(ch.col1));
279
-        if (ch.nChanges == 2)
280
-        {
281
-            a2 = alphaParameters_comp(D.nCol(), D.rowPtr(ch.row2), S.rowPtr(ch.row2),
282
-            AP.rowPtr(ch.row2), P.rowPtr(ch.col2));
283
-        }
284
-        return a1 + a2;
285
-    }
286
-    else if (ch.label == 'P' && ch.nChanges == 2 && ch.col1 == ch.col2)
287
-    {
288
-        return alphaParameters_comp(D.nRow(), D.colPtr(ch.col1), S.colPtr(ch.col1),
289
-            AP.colPtr(ch.col1), A.colPtr(ch.row1), A.colPtr(ch.row2));
290
-    }
291
-    else // single change, or two independent changes
206
+    for (unsigned i = 0; i < C.nRow(); ++i)
292 207
     {
293
-        AlphaParameters a1(0.f, 0.f), a2(0.f, 0.f);
294
-        a1 = alphaParameters_comp(D.nRow(), D.colPtr(ch.col1), S.colPtr(ch.col1),
295
-            AP.colPtr(ch.col1), A.colPtr(ch.row1));
296
-        if (ch.nChanges == 2)
208
+        for (unsigned j = 0; j < C.nCol(); ++j)
297 209
         {
298
-            a2 = alphaParameters_comp(D.nRow(), D.colPtr(ch.col2), S.colPtr(ch.col2),
299
-            AP.colPtr(ch.col2), A.colPtr(ch.row2));
210
+            float sum = 0.0;
211
+            for (unsigned k = 0; k < A.nCol(); ++k)
212
+            {
213
+                sum += A(i,k) * B(k,j);
214
+            }
215
+            C.set(i, j, sum);
300 216
         }
301
-        return a1 + a2;
302 217
     }
303 218
 }
304
-*/
305 219
\ No newline at end of file
220
+*/
... ...
@@ -26,6 +26,8 @@ namespace gaps
26 26
 {
27 27
 namespace algo
28 28
 {
29
+    const float epsilon = 1.0e-5f;
30
+
29 31
     // vector algorithms    
30 32
     unsigned whichMin(const Vector &vec);
31 33
     float sum(const Vector &vec);
... ...
@@ -50,7 +52,7 @@ namespace algo
50 52
         const GenericMatrix &meanMat, unsigned nUpdates);
51 53
 
52 54
     // specific matrix algorithms
53
-    bool isRowZero(const RowMatrix &mat, unsigned row);
55
+    bool isRowZero(const RowMatrix &mat, unsigned row); // TODO take pointer
54 56
     bool isColZero(const ColMatrix &mat, unsigned col);
55 57
     //void matrixMultiplication(TwoWayMatrix &C, const ColMatrix &A,
56 58
     //    const RowMatrix &B);
... ...
@@ -60,15 +62,19 @@ namespace algo
60 62
     float loglikelihood(const Matrix &D, const Matrix &S,
61 63
         const Matrix &AP);
62 64
 
63
-    // change in likelihood
64
-    //float deltaLL(const MatrixChange &ch, const TwoWayMatrix &D,
65
-    //    const TwoWayMatrix &S, const ColMatrix &A,
66
-    //    const RowMatrix &P, const TwoWayMatrix &AP);
65
+    AlphaParameters alphaParameters(unsigned size, const float *D,
66
+        const float *S, const float *AP, const float *mat);
67
+
68
+    AlphaParameters alphaParameters(unsigned size, const float *D,
69
+        const float *S, const float *AP, const float *mat1, const float *mat2);
70
+
71
+    float deltaLL(unsigned size, const float *D, const float *S,
72
+        const float *AP, const float *mat, float delta);
73
+
74
+    float deltaLL(unsigned size, const float *D, const float *S,
75
+        const float *AP, const float *mat1, float delta1, const float *mat2,
76
+        float delta2);
67 77
 
68
-    // alpha parameters used in exchange and gibbsMass calculation
69
-    //AlphaParameters alphaParameters(const MatrixChange &ch,
70
-    //    const TwoWayMatrix &D, const TwoWayMatrix &S, const ColMatrix &A,
71
-    //    const RowMatrix &P, const TwoWayMatrix &AP);*/
72 78
 } // namespace algo
73 79
 } // namespace gaps
74 80
 
... ...
@@ -5,11 +5,6 @@
5 5
 #include <stdint.h>
6 6
 #include <utility>
7 7
 
8
-unsigned AtomicDomain::size() const
9
-{
10
-    return mAtoms.size();
11
-}
12
-
13 8
 // O(1)
14 9
 Atom AtomicDomain::front() const
15 10
 {
... ...
@@ -19,8 +14,10 @@ Atom AtomicDomain::front() const
19 14
 // O(1)
20 15
 Atom AtomicDomain::randomAtom() const
21 16
 {
22
-    GAPS_ASSERT(!mAtoms.empty());
17
+    GAPS_ASSERT(mAtoms.size() > 0);
23 18
     uint64_t num = gaps::random::uniform64(0, mAtoms.size() - 1);
19
+    GAPS_ASSERT(num >= 0);
20
+    GAPS_ASSERT(num < mAtoms.size());
24 21
     return mAtoms[num];
25 22
 }
26 23
 
... ...
@@ -35,10 +32,57 @@ uint64_t AtomicDomain::randomFreePosition() const
35 32
     return pos;
36 33
 }
37 34
 
35
+// O(1)
36
+uint64_t AtomicDomain::size() const
37
+{
38
+    return mAtoms.size();
39
+}
40
+
41
+// O(1)
42
+Atom& AtomicDomain::left(const Atom &atom)
43
+{
44
+    GAPS_ASSERT(hasLeft(atom));
45
+    return mAtoms[atom.leftNdx - 1];
46
+}
47
+
48
+// O(1)
49
+Atom& AtomicDomain::right(const Atom &atom)
50
+{
51
+    GAPS_ASSERT(hasRight(atom));
52
+    return mAtoms[atom.rightNdx - 1];
53
+}
54
+
55
+// O(1)
56
+const Atom& AtomicDomain::left(const Atom &atom) const
57
+{
58
+    GAPS_ASSERT(hasLeft(atom));
59
+    return mAtoms[atom.leftNdx - 1];
60
+}
61
+
62
+// O(1)
63
+const Atom& AtomicDomain::right(const Atom &atom) const
64
+{
65
+    GAPS_ASSERT(hasRight(atom));
66
+    return mAtoms[atom.rightNdx - 1];
67
+}
68
+
69
+// O(1)
70
+bool AtomicDomain::hasLeft(const Atom &atom) const
71
+{
72
+    return atom.leftNdx != 0;
73
+}
74
+
75
+// O(1)
76
+bool AtomicDomain::hasRight(const Atom &atom) const
77
+{
78
+    return atom.rightNdx != 0;
79
+}
80
+
38 81
 // O(logN)
39 82
 void AtomicDomain::insert(uint64_t pos, float mass)
40 83
 {
41 84
     // insert position into map
85
+    GAPS_ASSERT(mAtomPositions.find(pos) == mAtomPositions.end());
42 86
     std::map<uint64_t, uint64_t>::iterator iter, iterLeft, iterRight;
43 87
     iter = mAtomPositions.insert(std::pair<uint64_t, uint64_t>(pos, mAtoms.size())).first;
44 88
     iterLeft = iter;
... ...
@@ -49,53 +93,106 @@ void AtomicDomain::insert(uint64_t pos, float mass)
49 93
     if (iter != mAtomPositions.begin())
50 94
     {
51 95
         --iterLeft;
52
-        atom.left = &(mAtoms[iterLeft->second]);
96
+        atom.leftNdx = iterLeft->second + 1;
53 97
     }
54 98
     if (++iter != mAtomPositions.end())
55 99
     {
56 100
         ++iterRight;
57
-        atom.right = &(mAtoms[iterRight->second]);
101
+        atom.rightNdx = iterRight->second + 1;
58 102
     } 
59 103
 
60 104
     // add atom to vector
105
+    GAPS_ASSERT(atom.rightNdx <= size());
106
+    GAPS_ASSERT(atom.leftNdx <= size());
61 107
     mAtoms.push_back(atom);
62 108
     mUsedPositions.insert(pos);
109
+    GAPS_ASSERT(mAtoms.size() == mAtomPositions.size());
110
+    GAPS_ASSERT(mAtoms.size() == mUsedPositions.size());
111
+}
112
+
113
+void AtomicDomain::swap(uint64_t i1, uint64_t i2)
114
+{
115
+    if (i1 != i2)
116
+    {
117
+        // store the atom in position 1
118
+        Atom temp1 = mAtoms[i1];
119
+        Atom temp2 = mAtoms[i2];
120
+
121
+        // switch atoms
122
+        mAtoms[i1] = temp2;
123
+        mAtoms[i2] = temp1;
124
+
125
+        // update neighbors
126
+        if (hasLeft(temp1))
127
+        {
128
+            GAPS_ASSERT(i2 + 1 <= size());
129
+            left(temp1).rightNdx = i2 + 1;
130
+        }
131
+        if (hasRight(temp1))
132
+        {
133
+            GAPS_ASSERT(i2 + 1 <= size());
134
+            right(temp1).leftNdx = i2 + 1;
135
+        }
136
+        if (hasLeft(temp2))
137
+        {
138
+            GAPS_ASSERT(i1 + 1 <= size());
139
+            left(temp2).rightNdx = i1 + 1;
140
+        }
141
+        if (hasRight(temp2))
142
+        {
143
+            GAPS_ASSERT(i1 + 1 <= size());
144
+            right(temp2).leftNdx = i1 + 1;
145
+        }
146
+
147
+        // update atom positions
148
+        mAtomPositions.erase(mAtoms[i1].pos);
149
+        mAtomPositions.erase(mAtoms[i2].pos);
150
+        mAtomPositions.insert(std::pair<uint64_t, uint64_t>(mAtoms[i1].pos,
151
+            i1));
152
+        mAtomPositions.insert(std::pair<uint64_t, uint64_t>(mAtoms[i2].pos,
153
+            i2));
154
+    }
63 155
 }
64 156
 
65 157
 // O(logN)
66
-void AtomicDomain::erase(uint64_t pos)
158
+// erase directly from position map
159
+// swap with last atom in vector, pop off the back
160
+void AtomicDomain::erase(uint64_t pos, bool display)
67 161
 {
68 162
     // get vector index of this atom and erase it
163
+    GAPS_ASSERT(mAtomPositions.find(pos) != mAtomPositions.end());
69 164
     uint64_t index = mAtomPositions.at(pos);
70
-    mAtomPositions.erase(pos);
71 165
 
72
-    // update key of object about to be moved (last one doesn't need to move)
73
-    if (index < mAtoms.size() - 1)
74
-    {
75
-        mAtomPositions.erase(mAtoms.back().pos);
76
-        mAtomPositions.insert(std::pair<uint64_t, uint64_t>(mAtoms.back().pos,
77
-            index));
78
-    }
166
+    // move atom to back
167
+    swap(index, mAtoms.size() - 1);
79 168
 
80
-    // update neighbors
81
-    if (mAtoms[index].left)
169
+    // connect neighbors of atom to be deleted
170
+    if (hasLeft(mAtoms.back()))
82 171
     {
83
-        mAtoms[index].left->right = mAtoms[index].right;
172
+        left(mAtoms.back()).rightNdx = mAtoms.back().rightNdx;
84 173
     }
85
-    if (mAtoms[index].right)
174
+    if (hasRight(mAtoms.back()))
86 175
     {
87
-        mAtoms[index].right->left = mAtoms[index].left;
176
+        right(mAtoms.back()).leftNdx = mAtoms.back().leftNdx;
88 177
     }
89 178
 
90 179
     // delete atom from vector in O(1)
91
-    mAtoms[index] = mAtoms.back();
180
+    GAPS_ASSERT(mAtomPositions.at(pos) == mAtoms.size() - 1);
181
+    mAtomPositions.erase(pos);
92 182
     mAtoms.pop_back();
93 183
     mUsedPositions.erase(pos);
184
+    GAPS_ASSERT(mAtoms.size() == mAtomPositions.size());
185
+    GAPS_ASSERT(mAtoms.size() == mUsedPositions.size());
94 186
 }
95 187
 
96 188
 // O(logN)
97 189
 void AtomicDomain::updateMass(uint64_t pos, float newMass)
98 190
 {
99
-    uint64_t index = mAtomPositions.at(pos);
100
-    mAtoms[index].mass = newMass;
191
+    GAPS_ASSERT(mAtomPositions.find(pos) != mAtomPositions.end());
192
+    mAtoms[mAtomPositions.at(pos)].mass = newMass;
101 193
 }
194
+
195
+void AtomicDomain::test(uint64_t pos) const
196
+{
197
+    GAPS_ASSERT(mAtomPositions.find(pos) != mAtomPositions.end());
198
+}
102 199
\ No newline at end of file
... ...
@@ -1,5 +1,5 @@
1
-#ifndef __GAPS_ATOMIC_DOMAIN_H__
2
-#define __GAPS_ATOMIC_DOMAIN_H__
1
+#ifndef __COGAPS_ATOMIC_DOMAIN_H__
2
+#define __COGAPS_ATOMIC_DOMAIN_H__
3 3
 
4 4
 #include "Archive.h"
5 5
 
... ...
@@ -9,16 +9,27 @@
9 9
 #include <vector>
10 10
 #include <map>
11 11
 
12
+class AtomicDomain;
13
+
14
+// Want the neighbors to be pointers, but can't use pointers since vector
15
+// is resized, shifted integers allow for 0 to be "null" and 1 to be the
16
+// first index. AtomicDomain is responsible for managing this correctly.
12 17
 struct Atom
13 18
 {
19
+//private:    
20
+public:
21
+
22
+    friend AtomicDomain;
23
+    uint64_t leftNdx; // shift these up 1, 0 == no neighbor
24
+    uint64_t rightNdx;
25
+
26
+public:    
27
+
14 28
     uint64_t pos;
15 29
     float mass;
16 30
 
17
-    Atom* left;
18
-    Atom* right;
19
-    
20 31
     Atom(uint64_t p, float m)
21
-        : pos(p), mass(m), left(NULL), right(NULL)
32
+        : pos(p), mass(m), leftNdx(0), rightNdx(0)
22 33
     {}
23 34
 
24 35
     bool operator==(const Atom &other) const
... ...
@@ -45,12 +56,21 @@ public:
45 56
     Atom front() const;
46 57
     Atom randomAtom() const;
47 58
     uint64_t randomFreePosition() const;
48
-    unsigned size() const;
59
+    uint64_t size() const;
60
+
61
+    Atom& left(const Atom &atom);
62
+    Atom& right(const Atom &atom);
63
+    const Atom& left(const Atom &atom) const;
64
+    const Atom& right(const Atom &atom) const;
65
+    bool hasLeft(const Atom &atom) const;
66
+    bool hasRight(const Atom &atom) const;
49 67
 
50 68
     // modify domain
51 69
     void insert(uint64_t pos, float mass);
52
-    void erase(uint64_t pos);
70
+    void swap(uint64_t i1, uint64_t i2);
71
+    void erase(uint64_t pos, bool display=false);
53 72
     void updateMass(uint64_t pos, float newMass);
73
+    void test(uint64_t pos) const;
54 74
 
55 75
     // serialization
56 76
     friend Archive& operator<<(Archive &ar, AtomicDomain &domain);
... ...
@@ -29,7 +29,7 @@ const std::string &cptFile, unsigned pumpThreshold, unsigned nPumpSamples)
29 29
 
30 30
     // create internal state from parameters and run from there
31 31
     GapsRunner runner(D, S, nFactor, nEquil, nEquilCool, nSample,
32
-        nOutputs, nSnapshots, alphaA, alphaP, maxGibbmassA, maxGibbmassP, seed,
32
+        nOutputs, nSnapshots, alphaA, alphaP, maxGibbmassA, maxGibbmassP, seedUsed,
33 33
         messages, singleCellRNASeq,  checkpointInterval, cptFile,
34 34
         whichMatrixFixed, FP);
35 35
     return runner.run();
... ...
@@ -1,5 +1,5 @@
1
-#ifndef __GAPS_ASSERT_H__
2
-#define __GAPS_ASSERT_H__
1
+#ifndef __COGAPS_GAPS_ASSERT_H__
2
+#define __COGAPS_GAPS_ASSERT_H__
3 3
 
4 4
 #ifdef GAPS_DEBUG
5 5
     #define GAPS_ASSERT(cond)                                           \
... ...
@@ -10,8 +10,19 @@
10 10
                 std::exit(0);                                           \
11 11
             }                                                           \
12 12
         } while(0)
13
+
14
+    #define GAPS_ASSERT_MSG(cond, msg)                                  \
15
+        do {                                                            \
16
+            if (!(cond))                                                \
17
+            {                                                           \
18
+                Rcpp::Rcout << "assert failed " << __FILE__ << " " <<   \
19
+                    __LINE__ << " , " << msg << '\n';                   \
20
+                std::exit(0);                                           \
21
+            }                                                           \
22
+        } while(0)
13 23
 #else
14 24
     #define GAPS_ASSERT(cond) ((void)sizeof(cond))
25
+    #define GAPS_ASSERT_MSG(cond, msg) ((void)sizeof(cond))
15 26
 #endif 
16 27
 
17 28
 #endif
18 29
\ No newline at end of file
... ...
@@ -40,6 +40,9 @@ mStatistics(D.nrow(), D.ncol(), nFactor)
40 40
    //    >> mPrintMessages >> mCurrentIter >> mPhase >> mSeed
41 41
    //    >> mCheckpointInterval >> mCheckpointFile >> mNumUpdatesA
42 42
    //    >> mNumUpdatesP >> mASampler >> mPSampler >> mStatistics;
43
+
44
+    mASampler.sync(mPSampler);
45
+    mPSampler.sync(mASampler);
43 46
 }
44 47
 
45 48
 // execute the steps of the algorithm, return list to R
... ...
@@ -95,6 +98,7 @@ void GapsRunner::runBurnPhase()
95 98
 {
96 99
     for (; mCurrentIter < mEquilIter; ++mCurrentIter)
97 100
     {
101
+        Rcpp::checkUserInterrupt();
98 102
         makeCheckpointIfNeeded();
99 103
         float temp = ((float)mCurrentIter + 2.f) / ((float)mEquilIter / 2.f);
100 104
         mASampler.setAnnealingTemp(std::min(1.f,temp));
... ...
@@ -109,6 +113,7 @@ void GapsRunner::runCoolPhase()
109 113
 {
110 114
     for (; mCurrentIter < mCoolIter; ++mCurrentIter)
111 115
     {
116
+        Rcpp::checkUserInterrupt();
112 117
         makeCheckpointIfNeeded();
113 118
         updateSampler();
114 119
     }
... ...
@@ -118,6 +123,7 @@ void GapsRunner::runSampPhase()
118 123
 {
119 124
     for (; mCurrentIter < mSampleIter; ++mCurrentIter)
120 125
     {
126
+        Rcpp::checkUserInterrupt();
121 127
         makeCheckpointIfNeeded();
122 128
         updateSampler();
123 129
         mStatistics.update(mASampler, mPSampler);
... ...
@@ -150,7 +156,7 @@ void GapsRunner::displayStatus(const std::string &type, unsigned nIterTotal)
150 156
 {
151 157
     if ((mCurrentIter + 1) % mNumOutputs == 0 && mPrintMessages)
152 158
     {
153
-        Rprintf("%s %d of %d, Atoms:%d(%d) Chi2 = %.2f\n", type.c_str(),
159
+        Rprintf("%s %d of %d, Atoms:%lu(%lu) Chi2 = %.2f\n", type.c_str(),
154 160
             mCurrentIter + 1, nIterTotal, mASampler.nAtoms(),
155 161
             mPSampler.nAtoms(), mASampler.chi2());
156 162
     }
... ...
@@ -69,13 +69,15 @@ private:
69 69
 
70 70
 public:
71 71
 
72
+    // construct from parameters
72 73
     GapsRunner(const Rcpp::NumericMatrix &D, const Rcpp::NumericMatrix &S,
73 74
         unsigned nFactor, unsigned nEquil, unsigned nCool, unsigned nSample,
74 75
         unsigned nOutputs, unsigned nSnapshots, float alphaA, float alphaP,
75 76
         float maxGibbsMassA, float maxGibbsMassP, uint32_t seed, bool messages,
76 77
         bool singleCellRNASeq, unsigned cptInterval, const std::string &cptFile,
77 78
         char whichMatrixFixed, const Rcpp::NumericMatrix &FP);
78
-
79
+    
80
+    // construct from checkpoint file
79 81
     GapsRunner(const Rcpp::NumericMatrix &D, const Rcpp::NumericMatrix &S,
80 82
         unsigned nFactor, unsigned nEquil, unsigned nSample,
81 83
         const std::string &cptFile);
... ...
@@ -15,16 +15,16 @@ const PatternGibbsSampler &PSampler)
15 15
     Vector normVec(mNumPatterns);
16 16
     for (unsigned j = 0; j < mNumPatterns; ++j)
17 17
     {
18
-        //normVec[j] = gaps::algo::sum(PSampler.mPMatrix.getRow(j));
19
-        //normVec[j] = normVec[j] == 0 ? 1.f : normVec[j];
18
+        normVec[j] = gaps::algo::sum(PSampler.mMatrix.getRow(j));
19
+        normVec[j] = normVec[j] == 0 ? 1.f : normVec[j];
20 20
 
21
-        //Vector quot(PSampler.mPMatrix.getRow(j) / normVec[j]);
22
-        //mPMeanMatrix.getRow(j) += quot;
23
-        //mPStdMatrix.getRow(j) += gaps::algo::elementSq(quot);
21
+        Vector quot(PSampler.mMatrix.getRow(j) / normVec[j]);
22
+        mPMeanMatrix.getRow(j) += quot;
23
+        mPStdMatrix.getRow(j) += gaps::algo::elementSq(quot);
24 24
 
25
-        //Vector prod(ASampler.mAMatrix.getCol(j) * normVec[j]);
26
-        //mAMeanMatrix.getCol(j) += prod;
27
-        //mAStdMatrix.getCol(j) += gaps::algo::elementSq(prod); 
25
+        Vector prod(ASampler.mMatrix.getCol(j) * normVec[j]);
26
+        mAMeanMatrix.getCol(j) += prod;
27
+        mAStdMatrix.getCol(j) += gaps::algo::elementSq(prod); 
28 28
     }
29 29
 }
30 30
 
... ...
@@ -36,7 +36,7 @@ bool AmplitudeGibbsSampler::canUseGibbs(unsigned r1, unsigned c1, unsigned r2, u
36 36
 
37 37
 void AmplitudeGibbsSampler::sync(PatternGibbsSampler &sampler)
38 38
 {
39
-    mOtherMatrix = &sampler.mMatrix;
39
+    mOtherMatrix = &(sampler.mMatrix);
40 40
     mAPMatrix = sampler.mAPMatrix;
41 41
 }
42 42
 
... ...
@@ -44,7 +44,50 @@ void AmplitudeGibbsSampler::updateAPMatrix(unsigned row, unsigned col, float del
44 44
 {
45 45
     for (unsigned j = 0; j < mAPMatrix.nCol(); ++j)
46 46
     {
47
-        mAPMatrix(row,j) += delta * (*mOtherMatrix)(j,col);
47
+        mAPMatrix(row,j) += delta * (*mOtherMatrix)(col,j);
48
+    }
49
+}
50
+
51
+AlphaParameters AmplitudeGibbsSampler::alphaParameters(unsigned row, unsigned col)
52
+{
53
+    return gaps::algo::alphaParameters(mDMatrix.nCol(), mDMatrix.rowPtr(row),
54
+        mSMatrix.rowPtr(row), mAPMatrix.rowPtr(row), mOtherMatrix->rowPtr(col));
55
+}
56
+
57
+AlphaParameters AmplitudeGibbsSampler::alphaParameters(unsigned r1, unsigned c1,
58
+unsigned r2, unsigned c2)
59
+{
60
+    if (r1 == r2)
61
+    {
62
+        return gaps::algo::alphaParameters(mDMatrix.nCol(), mDMatrix.rowPtr(r1),
63
+            mSMatrix.rowPtr(r1), mAPMatrix.rowPtr(r1), mOtherMatrix->rowPtr(c1),
64
+            mOtherMatrix->rowPtr(c2));
65
+    }
66
+    else
67
+    {
68
+        return alphaParameters(r1, c1) + alphaParameters(r2, c2);
69
+    }
70
+}
71
+
72
+float AmplitudeGibbsSampler::computeDeltaLL(unsigned row, unsigned col, float mass)
73
+{
74
+    return gaps::algo::deltaLL(mDMatrix.nCol(), mDMatrix.rowPtr(row),
75
+        mSMatrix.rowPtr(row), mAPMatrix.rowPtr(row), mOtherMatrix->rowPtr(col),
76
+        mass);
77
+}
78
+
79
+float AmplitudeGibbsSampler::computeDeltaLL(unsigned r1, unsigned c1, float m1,
80
+unsigned r2, unsigned c2, float m2)
81
+{
82
+    if (r1 == r2)
83
+    {
84
+        return gaps::algo::deltaLL(mDMatrix.nCol(), mDMatrix.rowPtr(r1),
85
+            mSMatrix.rowPtr(r1), mAPMatrix.rowPtr(r1), mOtherMatrix->rowPtr(c1),
86
+            m1, mOtherMatrix->rowPtr(c2), m2);
87
+    }
88
+    else
89
+    {
90
+        return computeDeltaLL(r1, c1, m1) + computeDeltaLL(r2, c2, m2);
48 91
     }
49 92
 }
50 93
 
... ...
@@ -84,7 +127,7 @@ bool PatternGibbsSampler::canUseGibbs(unsigned r1, unsigned c1, unsigned r2, uns
84 127
 
85 128
 void PatternGibbsSampler::sync(AmplitudeGibbsSampler &sampler)
86 129
 {
87
-    mOtherMatrix = &sampler.mMatrix;
130
+    mOtherMatrix = &(sampler.mMatrix);
88 131
     mAPMatrix = sampler.mAPMatrix;
89 132
 }
90 133
 
... ...
@@ -94,4 +137,47 @@ void PatternGibbsSampler::updateAPMatrix(unsigned row, unsigned col, float delta
94 137
     {
95 138
         mAPMatrix(i,col) += delta * (*mOtherMatrix)(i,row);
96 139
     }
140
+}
141
+
142
+AlphaParameters PatternGibbsSampler::alphaParameters(unsigned row, unsigned col)
143
+{
144
+    return gaps::algo::alphaParameters(mDMatrix.nRow(), mDMatrix.colPtr(col),
145
+        mSMatrix.colPtr(col), mAPMatrix.colPtr(col), mOtherMatrix->colPtr(row));
146
+}
147
+
148
+AlphaParameters PatternGibbsSampler::alphaParameters(unsigned r1, unsigned c1,
149
+unsigned r2, unsigned c2)
150
+{
151
+    if (c1 == c2)
152
+    {
153
+        return gaps::algo::alphaParameters(mDMatrix.nRow(), mDMatrix.colPtr(c1),
154
+            mSMatrix.colPtr(c1), mAPMatrix.colPtr(c1), mOtherMatrix->colPtr(r1),
155
+            mOtherMatrix->colPtr(r2));
156
+    }
157
+    else
158
+    {
159
+        return alphaParameters(r1, c1) + alphaParameters(r2, c2);
160
+    }
161
+}
162
+
163
+float PatternGibbsSampler::computeDeltaLL(unsigned row, unsigned col, float mass)
164
+{
165
+    return gaps::algo::deltaLL(mDMatrix.nRow(), mDMatrix.colPtr(col),
166
+        mSMatrix.colPtr(col), mAPMatrix.colPtr(col), mOtherMatrix->colPtr(row),
167
+        mass);
168
+}
169
+
170
+float PatternGibbsSampler::computeDeltaLL(unsigned r1, unsigned c1, float m1,
171
+unsigned r2, unsigned c2, float m2)
172
+{
173
+    if (c1 == c2)
174
+    {
175
+        return gaps::algo::deltaLL(mDMatrix.nRow(), mDMatrix.colPtr(c1),
176
+            mSMatrix.colPtr(c1), mAPMatrix.colPtr(c1), mOtherMatrix->colPtr(r1),
177
+            m1, mOtherMatrix->colPtr(r2), m2);
178
+    }
179
+    else
180
+    {
181
+        return computeDeltaLL(r1, c1, m1) + computeDeltaLL(r2, c2, m2);
182
+    }
97 183
 }
98 184
\ No newline at end of file
... ...
@@ -24,6 +24,7 @@ class GibbsSampler
24 24
 private:
25 25
 
26 26
     friend T; // prevent incorrect inheritance - only T can construct
27
+    friend GapsStatistics;
27 28
 
28 29
     GibbsSampler(const Rcpp::NumericMatrix &D, const Rcpp::NumericMatrix &S,
29 30
         unsigned nrow, unsigned ncol, float alpha);
... ...
@@ -50,11 +51,20 @@ protected:
50 51
     T* impl();
51 52
 
52 53
     void processProposal(const AtomicProposal &prop);
53
-    void birth(uint64_t pos);
54
-    void death(uint64_t pos, float mass);
55
-    void move(uint64_t src, float mass, uint64_t dest);
56
-    void exchange(uint64_t p1, float mass1, uint64_t p2, float mass2);
57
-    float gibbsMass(unsigned row, unsigned col);
54
+    void birth(uint64_t pos, unsigned row, unsigned col);
55
+    void death(uint64_t pos, float mass, unsigned row, unsigned col);
56
+    void move(uint64_t src, float mass, uint64_t dest, unsigned r1, unsigned c1,
57
+        unsigned r2, unsigned c2);
58
+    void exchange(uint64_t p1, float mass1, uint64_t p2, float mass2,
59
+        unsigned r1, unsigned c1, unsigned r2, unsigned c2);
60
+    float gibbsMass(unsigned row, unsigned col, float mass);
61
+    float gibbsMass(unsigned r1, unsigned c1, float m1, unsigned r2,
62
+        unsigned c2, float m2);
63
+    void addMass(uint64_t pos, float mass, unsigned row, unsigned col);
64
+    void removeMass(uint64_t pos, float mass, unsigned row, unsigned col);
65
+    void acceptExchange(uint64_t p1, float m1, float d1, uint64_t p2, float m2,
66
+        float d2, unsigned r1, unsigned c1, unsigned r2, unsigned c2);
67
+    bool updateAtomMass(uint64_t pos, float newMass);
58 68
 
59 69
 public:
60 70
 
... ...
@@ -62,11 +72,7 @@ public:
62 72
     void setAnnealingTemp(float temp);
63 73
     
64 74
     float chi2() const;
65
-    float nAtoms() const;
66
-
67
-    // serialization
68
-    //friend Archive& operator<<(Archive &ar, GibbsSampler &sampler);
69
-    //friend Archive& operator>>(Archive &ar, GibbsSampler &sampler);
75
+    uint64_t nAtoms() const;
70 76
 };
71 77
 
72 78
 class AmplitudeGibbsSampler : public GibbsSampler<AmplitudeGibbsSampler, ColMatrix, RowMatrix>
... ...
@@ -89,6 +95,14 @@ public:
89 95
         float maxGibbsMass=0.f);
90 96
 
91 97
     void sync(PatternGibbsSampler &sampler);
98
+
99
+    AlphaParameters alphaParameters(unsigned row, unsigned col);
100
+    AlphaParameters alphaParameters(unsigned r1, unsigned c1, unsigned r2,
101
+        unsigned c2);
102
+
103
+    float computeDeltaLL(unsigned row, unsigned col, float mass);
104
+    float computeDeltaLL(unsigned r1, unsigned c1, float m1, unsigned r2,
105
+        unsigned c2, float m2);
92 106
 };
93 107
 
94 108
 class PatternGibbsSampler : public GibbsSampler<PatternGibbsSampler, RowMatrix, ColMatrix>
... ...
@@ -111,6 +125,14 @@ public:
111 125
         float maxGibbsMass=0.f);
112 126
 
113 127
     void sync(AmplitudeGibbsSampler &sampler);
128
+
129
+    AlphaParameters alphaParameters(unsigned row, unsigned col);
130
+    AlphaParameters alphaParameters(unsigned r1, unsigned c1, unsigned r2,
131
+        unsigned c2);
132
+
133
+    float computeDeltaLL(unsigned row, unsigned col, float mass);
134
+    float computeDeltaLL(unsigned r1, unsigned c1, float m1, unsigned r2,
135
+        unsigned c2, float m2);
114 136
 };
115 137
 
116 138
 /******************* IMPLEMENTATION OF TEMPLATED FUNCTIONS *******************/
... ...
@@ -167,87 +189,237 @@ template <class T, class MatA, class MatB>
167 189
 void GibbsSampler<T, MatA, MatB>::processProposal(const AtomicProposal &prop)
168 190
 {
169 191
     GAPS_ASSERT(prop.type == 'B' || prop.type == 'D' || prop.type == 'M' || prop.type == 'E');
192
+    unsigned r1 = impl()->getRow(prop.pos1);
193
+    unsigned c1 = impl()->getCol(prop.pos1);
194
+    unsigned r2 = 0, c2 = 0;
170 195
     switch (prop.type)
171 196
     {
172
-        case 'B': birth(prop.pos1); break;
173
-        case 'D': death(prop.pos1, prop.mass1); break;
174
-        case 'M': move(prop.pos1, prop.mass1, prop.pos2); break;
175
-        case 'E': exchange(prop.pos1, prop.mass1, prop.pos2, prop.mass2); break;
197
+        case 'B':
198
+            birth(prop.pos1, r1, c1);
199
+            break;
200
+        case 'D':
201
+            death(prop.pos1, prop.mass1, r1, c1);
202
+            break;
203
+        case 'M':
204
+            r2 = impl()->getRow(prop.pos2);
205
+            c2 = impl()->getCol(prop.pos2);
206
+            move(prop.pos1, prop.mass1, prop.pos2, r1, c1, r2, c2);
207
+            break;
208
+        case 'E':
209
+            r2 = impl()->getRow(prop.pos2);
210
+            c2 = impl()->getCol(prop.pos2);
211
+            exchange(prop.pos1, prop.mass1, prop.pos2, prop.mass2, r1, c1, r2, c2);
212
+            break;
176 213
     }
177 214
 }
178 215
 
179 216
 template <class T, class MatA, class MatB>
180
-void GibbsSampler<T, MatA, MatB>::birth(uint64_t pos)
217
+void GibbsSampler<T, MatA, MatB>::addMass(uint64_t pos, float mass, unsigned row, unsigned col)
181 218
 {
182
-    unsigned row = impl()->getRow(pos);
183
-    unsigned col = impl()->getCol(pos);
184
-    float mass = impl()->canUseGibbs(row, col) ? gibbsMass(row, col)
185
-        : gaps::random::exponential(mLambda);
186
-
187 219
     mDomain.insert(pos, mass);
188 220
     mMatrix(row, col) += mass;
189
-    //impl()->updateAPMatrix(row, col, mass);
221
+    impl()->updateAPMatrix(row, col, mass);
222
+}
223
+template <class T, class MatA, class MatB>
224
+void GibbsSampler<T, MatA, MatB>::removeMass(uint64_t pos, float mass, unsigned row, unsigned col)
225
+{
226
+    mDomain.erase(pos);
227
+    mMatrix(row, col) += -mass;
228
+    impl()->updateAPMatrix(row, col, -mass);
190 229
 }
191 230
 
231
+// add an atom at pos, calculate mass either with an exponential distribution
232
+// or with the gibbs mass calculation
192 233
 template <class T, class MatA, class MatB>
193
-void GibbsSampler<T, MatA, MatB>::death(uint64_t pos, float mass)
234
+void GibbsSampler<T, MatA, MatB>::birth(uint64_t pos, unsigned row,
235
+unsigned col)
194 236
 {
195
-    mQueue.rejectDeath();
237
+    float mass = impl()->canUseGibbs(row, col) ? gibbsMass(row, col, 0.f)
238
+        : gaps::random::exponential(mLambda);
239
+    addMass(pos, mass, row, col);
196 240
 }
197 241
 
242
+// automatically accept death, attempt a rebirth at the same position, using
243
+// the original mass or the gibbs mass calculation
198 244
 template <class T, class MatA, class MatB>
199
-void GibbsSampler<T, MatA, MatB>::move(uint64_t src, float mass, uint64_t dest)
245
+void GibbsSampler<T, MatA, MatB>::death(uint64_t pos, float mass, unsigned row,
246
+unsigned col)
200 247
 {
201
-    unsigned r1 = impl()->getRow(src);
202
-    unsigned c1 = impl()->getCol(src);
203
-    unsigned r2 = impl()->getRow(dest);
204
-    unsigned c2 = impl()->getCol(dest);
205
-    if (r1 == r2 && c1 == c2)
248
+    removeMass(pos, mass, row, col);
249
+    mass = impl()->canUseGibbs(row, col) ? gibbsMass(row, col, mass) : mass;
250
+    float deltaLL = impl()->computeDeltaLL(row, col, mass);
251
+    if (deltaLL * mAnnealingTemp > std::log(gaps::random::uniform()))
206 252
     {
207
-        mDomain.erase(src);
208
-        mDomain.insert(dest, mass);
253
+        addMass(pos, mass, row, col);
254
+        mQueue.rejectDeath();
209 255
     }
210 256
     else
211 257
     {
212
-/*
213
-        if (deltaLL * mAnnealingTemp >= std::log(gaps::random::uniform()))
258
+        mQueue.acceptDeath();
259
+    }
260
+}
261
+
262
+// move mass from src to dest in the atomic domain
263
+template <class T, class MatA, class MatB>
264
+void GibbsSampler<T, MatA, MatB>::move(uint64_t src, float mass, uint64_t dest,
265
+unsigned r1, unsigned c1, unsigned r2, unsigned c2)
266
+{
267
+    if (r1 != r2 || c1 != c2) // automatically reject if change in same bin
268
+    {
269
+        float deltaLL = impl()->computeDeltaLL(r1, c1, -mass, r2, c2, mass);
270
+        if (deltaLL * mAnnealingTemp > std::log(gaps::random::uniform()))
214 271
         {
215
-            mDomain.deleteAtom(p1);
216
-            mDomain.addAtom(p2, mass);
217
-            mMatrix(r1, c1) += -mass;
218
-            mMatrix(r2, c2) += mass;
219
-            impl()->updateAPMatrix(r1, c1, -mass);
220
-            impl()->updateAPMatrix(r2, c2, mass);
272
+            removeMass(src, mass, r1, c1);
273
+            addMass(dest, mass, r2, c2);
221 274
         }
222
-*/
223 275
     }
224 276
 }
225 277
 
278
+// exchange some amount of mass between two positions, note it is possible
279
+// for one of the atoms to be deleted if it's mass becomes too small after
280
+// the exchange
226 281
 template <class T, class MatA, class MatB>
227
-void GibbsSampler<T, MatA, MatB>::exchange(uint64_t p1, float mass1, uint64_t p2, float mass2)
282
+void GibbsSampler<T, MatA, MatB>::exchange(uint64_t p1, float m1, uint64_t p2,
283
+float m2, unsigned r1, unsigned c1, unsigned r2, unsigned c2)
228 284
 {
285
+    if (r1 != r2 || c1 != c2) // automatically reject if change in same bin
286
+    {
287
+        if (m1 < m2)
288
+        {
289
+            if (impl()->canUseGibbs(r2, c2, r1, c1))
290
+            {
291
+                float gDelta = gibbsMass(r2, c2, m2, r1, c1, m1);
292
+                if (gDelta != 0.f)
293
+                {
294
+                    acceptExchange(p2, m2, gDelta, p1, m1, -gDelta, r2, c2, r1, c1);
295
+                    return;
296
+                }
297
+            }
298
+        }
299
+        else
300
+        {
301
+            if (impl()->canUseGibbs(r1, c1, r2, c2))
302
+            {
303
+                float gDelta = gibbsMass(r1, c1, m1, r2, c2, m2);
304
+                if (gDelta != 0.f)
305
+                {
306
+                    acceptExchange(p1, m1, gDelta, p2, m2, -gDelta, r1, c1, r2, c2);
307
+                    return;
308
+                }
309
+            }
310
+        }
311
+
312
+        // use metropolis hastings otherwise
313
+        float pUpper = gaps::random::p_gamma(m1 + m2, 2.f, 1.f / mLambda);
314
+        float newMass = gaps::random::inverseGammaSample(0.f, pUpper, 2.f, 1.f / mLambda);
315
+        float delta = m1 > m2 ? newMass - m1 : m2 - newMass; // change larger mass
316
+        float pOldMass = m1 + delta > m2 - delta ? m1 : m2;
317
+        float pNew = gaps::random::d_gamma(newMass, 2.f, 1.f / mLambda);
318
+        float pOld = gaps::random::d_gamma(pOldMass, 2.f, 1.f / mLambda);
319
+
320
+        if (pOld == 0.f && pNew != 0.f) // special case
321
+        {
322
+            acceptExchange(p1, m1, delta, p2, m2, -delta, r1, c1, r2, c2);
323
+            return;
324
+        }
325
+        else // normal case
326
+        {
327
+            float deltaLL = impl()->computeDeltaLL(r1, c1, delta, r2, c2, -delta);
328
+            float priorLL = (pOld == 0.f) ? 1.f : pOld / pNew;
329
+            float u = std::log(gaps::random::uniform() * priorLL);
330
+            if (u < deltaLL * mAnnealingTemp)
331
+            {
332
+                acceptExchange(p1, m1, delta, p2, m2, -delta, r1, c1, r2, c2);
333
+                return;
334
+            }
335
+        }
336
+    }
229 337
     mQueue.rejectDeath();
230 338
 }
231 339
 
340
+// helper function for acceptExchange, used to conditionally update the mass
341
+// at a single position, deleting the atom if the new mass is too small
342
+template <class T, class MatA, class MatB>
343
+bool GibbsSampler<T, MatA, MatB>::updateAtomMass(uint64_t pos, float newMass)
344
+{
345
+    if (newMass < gaps::algo::epsilon)
346
+    {
347
+        mDomain.erase(pos);
348
+        mQueue.acceptDeath();
349
+        return false;
350
+    }
351
+    else
352
+    {
353
+        mDomain.updateMass(pos, newMass);
354
+        return true;
355
+    }
356
+}
357
+
358
+// helper function for exchange step, updates the atomic domain, matrix, and
359
+// A*P matrix
232 360
 template <class T, class MatA, class MatB>
233
-float GibbsSampler<T, MatA, MatB>::gibbsMass(unsigned row, unsigned col)
361
+void GibbsSampler<T, MatA, MatB>::acceptExchange(uint64_t p1, float m1,
362
+float d1, uint64_t p2, float m2, float d2, unsigned r1, unsigned c1,
363
+unsigned r2, unsigned c2)
364
+{
365
+    bool b1 = updateAtomMass(p1, m1 + d1);
366
+    bool b2 = updateAtomMass(p2, m2 + d2);
367
+    GAPS_ASSERT(b1 || b2);
368
+
369
+    if (b1 && b2)
370
+    {
371
+        mQueue.rejectDeath();
372
+    }
373
+    mMatrix(r1, c1) += d1;
374
+    mMatrix(r2, c2) += d2;
375
+    impl()->updateAPMatrix(r1, c1, d1);
376
+    impl()->updateAPMatrix(r1, c2, d2);
377
+}
378
+
379
+template <class T, class MatA, class MatB>
380
+float GibbsSampler<T, MatA, MatB>::gibbsMass(unsigned row, unsigned col, float mass)
234 381
 {        
235
-    //AlphaParameters alpha = impl()->alphaParameters(row, col);
236
-    AlphaParameters alpha(10.f, 10.f);
237
-    alpha.s *= mAnnealingTemp / 2.f;
238
-    alpha.su *= mAnnealingTemp / 2.f;
239
-    float mean  = (2.f * alpha.su - mLambda) / (2.f * alpha.s);
240
-    float sd = 1.f / std::sqrt(2.f * alpha.s);
241
-
242
-    float plower = gaps::random::p_norm(0.f, mean, sd);
243
-
244
-    //float newMass = death ? mass : 0.f;
245
-    float newMass = 0.f;
246
-    if (plower < 1.f && alpha.s > 0.00001f)
382
+    AlphaParameters alpha = impl()->alphaParameters(row, col);
383
+    alpha.s *= mAnnealingTemp;
384
+    alpha.su *= mAnnealingTemp;
385
+
386
+    if (alpha.s > gaps::algo::epsilon)
247 387
     {
248
-        newMass = gaps::random::inverseNormSample(plower, 1.f, mean, sd);
388
+        float mean  = (alpha.su - mLambda) / alpha.s;
389
+        float sd = 1.f / std::sqrt(alpha.s);
390
+        float pLower = gaps::random::p_norm(0.f, mean, sd);
391
+    
392
+        if (pLower < 1.f)
393
+        {
394
+            mass = gaps::random::inverseNormSample(pLower, 1.f, mean, sd);
395
+        }
249 396
     }
250
-    return std::max(0.f, std::min(newMass, mMaxGibbsMass));
397
+    return std::min(std::max(gaps::algo::epsilon, mass), mMaxGibbsMass);
398
+}
399
+
400
+template <class T, class MatA, class MatB>
401
+float GibbsSampler<T, MatA, MatB>::gibbsMass(unsigned r1, unsigned c1, float m1,
402
+unsigned r2, unsigned c2, float m2)
403
+{
404
+    AlphaParameters alpha = impl()->alphaParameters(r1, c1, r2, c2);
405
+    alpha.s *= mAnnealingTemp;
406
+    alpha.su *= mAnnealingTemp;
407
+
408
+    if (alpha.s > gaps::algo::epsilon)
409
+    {
410
+        float mean = alpha.su / alpha.s;
411
+        float sd = 1.f / std::sqrt(alpha.s);
412
+        float pLower = gaps::random::p_norm(-m1, mean, sd);
413
+        float pUpper = gaps::random::p_norm(m2, mean, sd);
414
+
415
+        if (pLower <  0.95f && pUpper > 0.05f)
416
+        {
417
+            GAPS_ASSERT_MSG(pLower < pUpper, m1 << " , " << m2 << " , " << mean << " , " << sd);
418
+            float delta = gaps::random::inverseNormSample(pLower, pUpper, mean, sd);
419
+            return std::min(std::max(-m1, delta), m2); // conserve mass
420
+        }
421
+    }
422
+    return 0.f;
251 423
 }
252 424
 
253 425
 template <class T, class MatA, class MatB>
... ...
@@ -263,132 +435,9 @@ float GibbsSampler<T, MatA, MatB>::chi2() const
263 435
 }
264 436
   
265 437
 template <class T, class MatA, class MatB>
266
-float GibbsSampler<T, MatA, MatB>::nAtoms() const
438
+uint64_t GibbsSampler<T, MatA, MatB>::nAtoms() const
267 439
 {   
268 440
     return mDomain.size();
269 441
 }
270 442
 
271
-
272
-//
273
-
274
-//
275
-
276
-//
277
-////  template <class T, class MatA, class MatB>
278
-////  void GibbsSampler<T, MatA, MatB>::processProposal(const AtomicProposal &prop)
279
-////  {
280
-////      GAPS_ASSERT(prop.type == 'B' || prop.type == 'D' || prop.type == 'M' || prop.type == 'E');
281
-////      switch (prop.type)
282
-////      {
283
-////          case 'B': birth(prop.pos1); break;
284
-////          //case 'D': death(prop.pos1, prop.mass1); break;
285
-////          //case 'M': move(prop.pos1, prop.mass1, prop.pos2); break;
286
-////          //case 'E': exchange(prop.pos1, prop.mass1, prop.pos2, prop.mass2); break;
287
-////      }
288
-////  }
289
-////  
290
-//
291
-////  
292
-////  template <class T, class MatA, class MatB>
293
-////  void GibbsSampler<T, MatA, MatB>::death(uint64_t pos, float mass)
294
-////  {
295
-////      /*unsigned row = impl()->getRow(pos);
296
-////      unsigned col = impl()->getCol(pos);
297
-////      mMatrix(row, col) += -mass;
298
-////      impl()->updateAPMatrix(row, col, -mass);
299
-////      mDomain.erase(pos);
300
-////      mQueue.acceptDeath();
301
-////  
302
-////      float newMass = impl()->canUseGibbs(row, col) ? gibbsMass(row, col)
303
-////      : mass;
304
-//      float deltaLL = impl()->computeDeltaLL(row, col, newMass);
305
-//  
306
-//      if (deltaLL * mAnnealingTemp >= std::log(gaps::random::uniform()))
307
-//      {
308
-//      float delta = mDomain.updateAtomMass(pos, newMass - mass);
309
-//      mMatrix(row, col) += delta;
310
-//      impl()->updateAPMatrix(row, col, delta);
311
-//      if (mDomain.at(pos))
312
-//      {
313
-//          mQueue.rejectDeath();
314
-//      }
315
-//      }
316
-//      else
317
-//      {
318
-//      mDomain.deleteAtom(pos);
319
-//      mQueue.maxAtoms--;
320
-//      }*/
321
-//  }
322
-//  
323
-//  template <class T, class MatA, class MatB>
324
-//  void GibbsSampler<T, MatA, MatB>::move(uint64_t src, float mass, uint64_t dest)
325
-//  {
326
-//      /*
327
-//      if (r1 == r2 && c1 == c2)
328
-//      {
329
-//      mDomain.deleteAtom(p1);
330
-//      mDomain.addAtom(p2, mass);
331
-//      }
332
-//      else
333
-//      {
334
-//      if (deltaLL * mAnnealingTemp >= std::log(gaps::random::uniform()))
335
-//      {
336
-//          mDomain.deleteAtom(p1);
337
-//          mDomain.addAtom(p2, mass);
338
-//          mMatrix(r1, c1) += -mass;
339
-//          mMatrix(r2, c2) += mass;
340
-//          impl()->updateAPMatrix(r1, c1, -mass);
341
-//          impl()->updateAPMatrix(r2, c2, mass);
342
-//      }
343
-//      }
344
-//      */
345
-//  }
346
-//  
347
-//  template <class T, class MatA, class MatB>
348
-//  void GibbsSampler<T, MatA, MatB>::exchange(uint64_t p1, float mass1, uint64_t p2, float mass2)
349
-//  {
350
-//      /*
351
-//      float newMass = gaps::random::inverseGammaSample(0.f, mass1 + mass2, 2.f, 1.f / mLambda);
352
-//      float delta1 = mass1 > mass2 ? newMass - mass1 : mass2 - newMass;
353
-//      float delta2 = mass2 > mass1 ? newMass - mass2 : mass1 - newMass;
354
-//      unsigned r1 = impl()->getRow(p1);
355
-//      unsigned c1 = impl()->getCol(p1);
356
-//      unsigned r2 = impl()->getRow(p2);
357
-//      unsigned c2 = impl()->getCol(p2);
358
-//  
359
-//      if (r1 == r2 && c1 == c2)
360
-//      {
361
-//      mDomain.updateAtomMass(p1, delta1);
362
-//      mDomain.updateAtomMass(p2, delta2);
363
-//      }
364
-//      else
365
-//      {
366
-//      // all the exchange code
367
-//      }
368
-//      */
369
-//      //mQueue.rejectDeath();
370
-//  }
371
-//  
372
-//  // don't return mass less than epsilon
373
-//  template <class T, class MatA, class MatB>
374
-//  float GibbsSampler<T, MatA, MatB>::gibbsMass(unsigned row, unsigned col)
375
-//  {        
376
-//  /*
377
-//      AlphaParameters alpha = impl()->alphaParameters(row, col);
378
-//      alpha.s *= mAnnealingTemp / 2.f;
379
-//      alpha.su *= mAnnealingTemp / 2.f;
380
-//      float mean  = (2.f * alpha.su - mLambda) / (2.f * alpha.s);
381
-//      float sd = 1.f / std::sqrt(2.f * alpha.s);
382
-//  
383
-//      float plower = gaps::random::p_norm(0.f, mean, sd);
384
-//  
385
-//      float newMass = death ? mass : 0.f;
386
-//      if (plower < 1.f && alpha.s > 0.00001f)
387
-//      {
388
-//      newMass = gaps::random::inverseNormSample(plower, 1.f, mean, sd);
389
-//      }
390
-//      return std::max(0.f, std::min(newMax, mMaxGibbsMass));
391
-//  */
392
-//  }
393
-
394 443
 #endif
395 444
\ No newline at end of file
... ...
@@ -1,4 +1,4 @@
1
-PKG_CPPFLAGS = -DGAPS_DEBUG -DBOOST_MATH_PROMOTE_DOUBLE_POLICY=0 -DSIMD -msse4.1
1
+PKG_CPPFLAGS = -DGAPS_DEBUG -DBOOST_MATH_PROMOTE_DOUBLE_POLICY=0
2 2
 OBJECTS =   Algorithms.o \
3 3
             AtomicDomain.o \
4 4
             Cogaps.o \
... ...
@@ -2,9 +2,6 @@
2 2
 #include "ProposalQueue.h"
3 3
 #include "Random.h"
4 4
 
5
-const uint64_t atomicEnd = std::numeric_limits<uint64_t>::max();
6
-const double atomicSize = static_cast<double>(atomicEnd);
7
-
8 5
 void ProposalQueue::populate(const AtomicDomain &domain, unsigned limit)
9 6
 {
10 7
     unsigned nIter = 0;
... ...
@@ -59,9 +56,16 @@ void ProposalQueue::rejectDeath()
59 56
     mMinAtoms++;
60 57
 }
61 58
 
59
+void ProposalQueue::rejectBirth()
60
+{
61
+    mMinAtoms--;
62
+    mMaxAtoms--;
63
+}
64
+
62 65
 float ProposalQueue::deleteProb(unsigned nAtoms) const
63 66
 {
64
-    double term1 = (atomicSize - static_cast<double>(nAtoms)) / atomicSize;
67
+    double size = static_cast<double>(mDomainSize);
68
+    double term1 = (size - static_cast<double>(nAtoms)) / size;
65 69
     double term2 = mAlpha * static_cast<double>(mNumBins) * term1;
66 70
     return static_cast<double>(nAtoms) / (static_cast<double>(nAtoms) + term2);
67 71
 }
... ...
@@ -138,8 +142,8 @@ bool ProposalQueue::death(const AtomicDomain &domain)
138 142
 bool ProposalQueue::move(const AtomicDomain &domain)
139 143
 {
140 144
     Atom a = domain.randomAtom();
141
-    uint64_t lbound = a.left ? a.left->pos : 0;
142
-    uint64_t rbound = a.right ? a.right->pos : atomicEnd - 1;
145
+    uint64_t lbound = domain.hasLeft(a) ? domain.left(a).pos : 0;
146
+    uint64_t rbound = domain.hasRight(a) ? domain.right(a).pos : mDomainSize;
143 147
 
144 148
     //for (unsigned i = 0; i < mUsedPositions.size(); ++i)
145 149
     //{
... ...
@@ -149,7 +153,7 @@ bool ProposalQueue::move(const AtomicDomain &domain)
149 153
     //    }
150 154
     //}
151 155
 
152
-    uint64_t newLocation = gaps::random::uniform64(lbound, rbound);
156
+    uint64_t newLocation = gaps::random::uniform64(lbound + 1, rbound - 1);
153 157
     //if (isInVector(mUsedIndices, a.pos / mDimensionSize) || isInVector(mUsedIndices, newLocation / mDimensionSize))
154 158
     //{
155 159
     //    return false; // matrix conflict - can't compute deltaLL
... ...
@@ -166,7 +170,9 @@ bool ProposalQueue::move(const AtomicDomain &domain)
166 170
 bool ProposalQueue::exchange(const AtomicDomain &domain)
167 171
 {
168 172
     Atom a1 = domain.randomAtom();
169
-
173
+    domain.test(a1.pos);
174
+    GAPS_ASSERT(a1.rightNdx <= domain.size());
175
+    GAPS_ASSERT(a1.leftNdx <= domain.size());
170 176
     //if (a1.right) // has neighbor
171 177
     //{
172 178
     //    for (unsigned i = 0; i < mUsedPositions.size(); ++i)
... ...
@@ -188,7 +194,13 @@ bool ProposalQueue::exchange(const AtomicDomain &domain)
188 194
     //    }
189 195
     //}
190 196
 
191
-    Atom a2 = a1.right ? *a1.right : domain.front();
197
+    Atom a2 = domain.hasRight(a1) ? domain.right(a1) : domain.front();
198
+    domain.test(a1.pos);
199
+    GAPS_ASSERT(a1.rightNdx <= domain.size());
200
+    GAPS_ASSERT(a1.leftNdx <= domain.size());
201
+    GAPS_ASSERT(a2.rightNdx <= domain.size());
202
+    GAPS_ASSERT(a2.leftNdx <= domain.size());
203
+    domain.test(a2.pos);
192 204
     //if (isInVector(mUsedIndices, a1.pos / mDimensionSize) || isInVector(mUsedIndices, a2.pos / mDimensionSize))
193 205
     //{
194 206
     //    return false; // matrix conflict - can't compute gibbs mass or deltaLL
... ...
@@ -70,6 +70,7 @@ public:
70 70
     // update min/max atoms
71 71
     void acceptDeath();
72 72
     void rejectDeath();
73
+    void rejectBirth();
73 74
 
74 75
     // serialization
75 76
     friend Archive& operator<<(Archive &ar, ProposalQueue &queue);
... ...
@@ -1,6 +1,7 @@
1 1
 // [[Rcpp::depends(BH)]]
2 2
 
3 3
 #include "Random.h"
4
+#include "GapsAssert.h"
4 5
 
5 6
 #include <boost/random/uniform_01.hpp>
6 7
 #include <boost/random/uniform_real_distribution.hpp>
... ...
@@ -22,7 +23,6 @@ typedef boost::random::mt19937 RNGType;
22 23
 //typedef boost::random::mt11213b RNGType; // should be faster
23 24
 
24 25
 static RNGType rng;
25
-static boost::random::uniform_01<RNGType&> u01_dist(rng);
26 26
 
27 27
 void gaps::random::save(Archive &ar)
28 28
 {
... ...
@@ -60,6 +60,7 @@ float gaps::random::exponential(float lambda)
60 60
 // open interval
61 61
 float gaps::random::uniform()
62 62
 {
63
+    boost::random::uniform_01<RNGType&> u01_dist(rng);
63 64
     return u01_dist();
64 65
 }
65 66
 
... ...
@@ -141,6 +142,8 @@ float gaps::random::p_norm(float p, float mean, float sd)
141 142
 
142 143
 float gaps::random::inverseNormSample(float a, float b, float mean, float sd)
143 144
 {
145
+    //GAPS_ASSERT(a != b);
146
+    GAPS_ASSERT(!((a == 0.f && b == 0.f) || (a == 1.f) && (b == 1.f)));
144 147
     float u = gaps::random::uniform(a, b);
145 148
     while (u == 0.f || u == 1.f)
146 149
     {
... ...
@@ -151,6 +154,7 @@ float gaps::random::inverseNormSample(float a, float b, float mean, float sd)
151 154
 
152 155
 float gaps::random::inverseGammaSample(float a, float b, float mean, float sd)
153 156
 {
157
+    GAPS_ASSERT(a != b);
154 158
     float u = gaps::random::uniform(a, b);
155 159
     while (u == 0.f || u == 1.f)
156 160
     {