1. CoGAPS
  2. Commit 094c9c46c1110d9c30d1f2b8240ba8d7807be8d3
Browse code

separated gibbs samplers

sherman5 authored on 07/03/2018 22:59:17
Showing14 changed files
src/Algorithms.cpp
History View file @094c9c46c
... ... 
@@ -48,7 +48,7 @@ unsigned gaps::algo::whichMin(const Vector &vec)
48 48 
 {
49 49 
     float min = vec[0];
50 50 
     unsigned minNdx = 0;
51 
-    for (unsigned i = 0; i < vec.size(); ++i)
51 
+    for (unsigned i = 1; i < vec.size(); ++i)
52 52 
     {
53 53 
         if (vec[i] < min)
54 54 
         {
... ... 
@@ -94,6 +94,7 @@ bool gaps::algo::isColZero(const ColMatrix &mat, unsigned col)
94 94 
 {
95 95 
     return gaps::algo::sum(mat.getCol(col)) == 0;
96 96 
 }
97 
+
97 98 
 /*
98 99 
 // horribly slow, don't call often
99 100 
 void gaps::algo::matrixMultiplication(TwoWayMatrix &C, const ColMatrix &A,
... ... 
@@ -113,20 +114,6 @@ const RowMatrix &B)
113 114 
     }
114 115 
 }
115 116
116 
-float gaps::algo::loglikelihood(const TwoWayMatrix &D,
117 
-const TwoWayMatrix &S, const TwoWayMatrix &AP)
118 
-{
119 
-    float chi2 = 0.f;
120 
-    for (unsigned i = 0; i < D.nRow(); ++i)
121 
-    {
122 
-        for (unsigned j = 0; j < D.nCol(); ++j)
123 
-        {
124 
-            chi2 += GAPS_SQ(D(i,j) - AP(i,j)) / GAPS_SQ(S(i,j));
125 
-        }
126 
-    }
127 
-    return chi2 / 2.f;
128 
-}
129 
-
130 117 
 static float deltaLL_comp(unsigned size, const float *D, const float *S,
131 118 
 const float *AP, const float *other, float delta)
132 119 
 {
src/Algorithms.h
History View file @094c9c46c
... ... 
@@ -52,22 +52,23 @@ namespace algo
52 52 
     // specific matrix algorithms
53 53 
     bool isRowZero(const RowMatrix &mat, unsigned row);
54 54 
     bool isColZero(const ColMatrix &mat, unsigned col);
55 
-    /*void matrixMultiplication(TwoWayMatrix &C, const ColMatrix &A,
56 
-        const RowMatrix &B);
55 
+    //void matrixMultiplication(TwoWayMatrix &C, const ColMatrix &A,
56 
+    //    const RowMatrix &B);
57 57
58 58 
     // chiSq / 2
59 
-    float loglikelihood(const TwoWayMatrix &D, const TwoWayMatrix &S,
60 
-        const TwoWayMatrix &AP);
59 
+    template <class Matrix>
60 
+    float loglikelihood(const Matrix &D, const Matrix &S,
61 
+        const Matrix &AP);
61 62
62 63 
     // change in likelihood
63 
-    float deltaLL(const MatrixChange &ch, const TwoWayMatrix &D,
64 
-        const TwoWayMatrix &S, const ColMatrix &A,
65 
-        const RowMatrix &P, const TwoWayMatrix &AP);
64 
+    //float deltaLL(const MatrixChange &ch, const TwoWayMatrix &D,
65 
+    //    const TwoWayMatrix &S, const ColMatrix &A,
66 
+    //    const RowMatrix &P, const TwoWayMatrix &AP);
66 67
67 68 
     // alpha parameters used in exchange and gibbsMass calculation
68 
-    AlphaParameters alphaParameters(const MatrixChange &ch,
69 
-        const TwoWayMatrix &D, const TwoWayMatrix &S, const ColMatrix &A,
70 
-        const RowMatrix &P, const TwoWayMatrix &AP);*/
69 
+    //AlphaParameters alphaParameters(const MatrixChange &ch,
70 
+    //    const TwoWayMatrix &D, const TwoWayMatrix &S, const ColMatrix &A,
71 
+    //    const RowMatrix &P, const TwoWayMatrix &AP);*/
71 72 
 } // namespace algo
72 73 
 } // namespace gaps
73 74
... ... 
@@ -127,4 +128,19 @@ const GenericMatrix &meanMat, unsigned nUpdates)
127 128 
     return retMat;
128 129 
 }
129 130
131 
+template <class Matrix>
132 
+float gaps::algo::loglikelihood(const Matrix &D, const Matrix &S,
133 
+const Matrix &AP)
134 
+{
135 
+    float chi2 = 0.f;
136 
+    for (unsigned i = 0; i < D.nRow(); ++i)
137 
+    {
138 
+        for (unsigned j = 0; j < D.nCol(); ++j)
139 
+        {
140 
+            chi2 += GAPS_SQ(D(i,j) - AP(i,j)) / GAPS_SQ(S(i,j));
141 
+        }
142 
+    }
143 
+    return chi2 / 2.f;
144 
+}
145 
+
130 146 
 #endif
131 147 
\ No newline at end of file
src/AtomicDomain.cpp
History View file @094c9c46c
... ... 
@@ -1,9 +1,15 @@
1 
+#include "GapsAssert.h"
1 2 
 #include "AtomicDomain.h"
2 3 
 #include "Random.h"
3 4
4 5 
 #include <stdint.h>
5 6 
 #include <utility>
6 7
8 
+unsigned AtomicDomain::size() const
9 
+{
10 
+    return mAtoms.size();
11 
+}
12 
+
7 13 
 // O(1)
8 14 
 Atom AtomicDomain::front() const
9 15 
 {
... ... 
@@ -13,18 +19,19 @@ Atom AtomicDomain::front() const
13 19 
 // O(1)
14 20 
 Atom AtomicDomain::randomAtom() const
15 21 
 {
22 
+    GAPS_ASSERT(!mAtoms.empty());
16 23 
     uint64_t num = gaps::random::uniform64(0, mAtoms.size() - 1);
17 24 
     return mAtoms[num];
18 25 
 }
19 26
20 
-// O(logN) - keep hash of positions to fix, need this O(1)
27 
+// Average Case O(1)
21 28 
 uint64_t AtomicDomain::randomFreePosition() const
22 29 
 {
23 30 
     uint64_t pos = 0;
24 31 
     do
25 32 
     {
26 33 
         pos = gaps::random::uniform64();
27 
-    } while (mAtomPositions.count(pos) > 0); // count is O(logN)
34 
+    } while (mUsedPositions.count(pos) > 0); // hash map => count is O(l)
28 35 
     return pos;
29 36 
 }
30 37
... ... 
@@ -44,7 +51,7 @@ void AtomicDomain::insert(uint64_t pos, float mass)
44 51 
         --iterLeft;
45 52 
         atom.left = &(mAtoms[iterLeft->second]);
46 53 
     }
47 
-    if (iter != mAtomPositions.end())
54 
+    if (++iter != mAtomPositions.end())
48 55 
     {
49 56 
         ++iterRight;
50 57 
         atom.right = &(mAtoms[iterRight->second]);
... ... 
@@ -52,6 +59,7 @@ void AtomicDomain::insert(uint64_t pos, float mass)
52 59
53 60 
     // add atom to vector
54 61 
     mAtoms.push_back(atom);
62 
+    mUsedPositions.insert(pos);
55 63 
 }
56 64
57 65 
 // O(logN)
... ... 
@@ -82,6 +90,7 @@ void AtomicDomain::erase(uint64_t pos)
82 90 
     // delete atom from vector in O(1)
83 91 
     mAtoms[index] = mAtoms.back();
84 92 
     mAtoms.pop_back();
93 
+    mUsedPositions.erase(pos);
85 94 
 }
86 95
87 96 
 // O(logN)
src/AtomicDomain.h
History View file @094c9c46c
... ... 
@@ -3,6 +3,7 @@
3 3
4 4 
 #include "Archive.h"
5 5
6 
+#include <boost/unordered_set.hpp>
6 7 
 #include <stdint.h>
7 8 
 #include <cstddef>
8 9 
 #include <vector>
... ... 
@@ -35,11 +36,16 @@ private:
35 36 
     std::vector<Atom> mAtoms;
36 37 
     std::map<uint64_t, uint64_t> mAtomPositions;
37 38
39 
+    // TODO google_dense_set - first profile and benchmark
40 
+    boost::unordered_set<uint64_t> mUsedPositions;
41 
+
38 42 
 public:
39 43
44 
+    // access atoms
40 45 
     Atom front() const;
41 46 
     Atom randomAtom() const;
42 47 
     uint64_t randomFreePosition() const;
48 
+    unsigned size() const;
43 49
44 50 
     // modify domain
45 51 
     void insert(uint64_t pos, float mass);
src/Cogaps.cpp
History View file @094c9c46c
... ... 
@@ -1,3 +1,4 @@
1 
+#include "GapsAssert.h"
1 2 
 #include "GibbsSampler.h"
2 3 
 #include "Matrix.h"
3 4 
 #include "Archive.h"
... ... 
@@ -52,12 +53,12 @@ static void createCheckpoint(GapsInternalState &state)
52 53 
 static void updateSampler(GapsInternalState &state)
53 54 
 {
54 55 
     state.nUpdatesA += state.nIterA;
55 
-    state.ASampler.update(state.nIterA);
56 
-    state.PSampler.syncAP(state.ASampler.APMatrix());
56 
+    update(state.ASampler, state.nIterA);
57 
+    state.PSampler.sync(state.ASampler);
57 58
58 59 
     state.nUpdatesP += state.nIterP;
59 
-    state.PSampler.update(state.nIterP);
60 
-    state.ASampler.syncAP(state.PSampler.APMatrix());
60 
+    update(state.PSampler, state.nIterP);
61 
+    state.ASampler.sync(state.PSampler);
61 62 
 }
62 63
63 64 
 static void makeCheckpointIfNeeded(GapsInternalState &state)
... ... 
@@ -183,6 +184,10 @@ static Rcpp::List runCogaps(GapsInternalState &state)
183 184 
         //Rcpp::Named("Asd") = state.sampler.AStdRMatrix(),
184 185 
         //Rcpp::Named("Pmean") = state.sampler.PMeanRMatrix(),
185 186 
         //Rcpp::Named("Psd") = state.sampler.PStdRMatrix(),
187 
+        Rcpp::Named("Amean") = Rcpp::NumericMatrix(193,9),
188 
+        Rcpp::Named("Asd") = Rcpp::NumericMatrix(193,9),
189 
+        Rcpp::Named("Pmean") = Rcpp::NumericMatrix(9,193),
190 
+        Rcpp::Named("Psd") = Rcpp::NumericMatrix(9,193),
186 191 
         Rcpp::Named("ASnapshots") = Rcpp::wrap(state.snapshotsA),
187 192 
         Rcpp::Named("PSnapshots") = Rcpp::wrap(state.snapshotsP),
188 193 
         Rcpp::Named("atomsAEquil") = state.nAtomsAEquil.rVec(),
src/GapsAssert.h
History View file @094c9c46c
189 194 
new file mode 100644
... ... 
@@ -0,0 +1,17 @@
1 
+#ifndef __GAPS_ASSERT_H__
2 
+#define __GAPS_ASSERT_H__
3 
+
4 
+#ifdef GAPS_DEBUG
5 
+    #define GAPS_ASSERT(cond)                                           \
6 
+        do {                                                            \
7 
+            if (!(cond))                                                \
8 
+            {                                                           \
9 
+                Rprintf("assert failed %s %d\n", __FILE__, __LINE__);   \
10 
+                std::exit(0);                                           \
11 
+            }                                                           \
12 
+        } while(0)
13 
+#else
14 
+    #define GAPS_ASSERT(cond) ((void)sizeof(cond))
15 
+#endif
16 
+
17 
+#endif
0 18 
\ No newline at end of file
src/GapsStatistics.cpp
History View file @094c9c46c
1 19 
new file mode 100644
... ... 
@@ -0,0 +1,38 @@
1 
+#include "GapsStatistics.h"
2 
+
3 
+GapsStatistics::GapsStatistics(unsigned nRow, unsigned nCol, unsigned nFactor)
4 
+    : mAMeanMatrix(nRow, nFactor), mAStdMatrix(nRow, nFactor),
5 
+        mPMeanMatrix(nFactor, nCol), mPStdMatrix(nFactor, nCol), mStatUpdates(0)
6 
+{}
7 
+
8 
+void GapsStatistics::update(const AmplitudeGibbsSampler &ASampler,
9 
+const PatternGibbsSampler &PSampler)
10 
+{
11 
+    mStatUpdates++;
12 
+
13 
+    mAMeanMatrix += ASampler.mAMatrix;
14 
+    mPMeanMatrix += PSampler.mPMatrix;
15 
+
16 
+}
17 
+
18 
+Rcpp::NumericMatrix GapsStatistics::AMean() const
19 
+{
20 
+    return (mAMeanMatrix / mStatUpdates).rMatrix();
21 
+}
22 
+
23 
+Rcpp::NumericMatrix GapsStatistics::AStd() const
24 
+{
25 
+    return gaps::algo::computeStdDev(mAStdMatrix, mAMeanMatrix,
26 
+        mStatUpdates).rMatrix();
27 
+}
28 
+
29 
+Rcpp::NumericMatrix GapsStatistics::PMean() const
30 
+{
31 
+    return (mPMeanMatrix / mStatUpdates).rMatrix();
32 
+}
33 
+
34 
+Rcpp::NumericMatrix GapsStatistics::PStd() const
35 
+{
36 
+    return gaps::algo::computeStdDev(mPStdMatrix, mPMeanMatrix,
37 
+        mStatUpdates).rMatrix();
38 
+}
src/GapsStatistics.h
History View file @094c9c46c
0 39 
new file mode 100644
... ... 
@@ -0,0 +1,31 @@
1 
+#ifndef __GAPS_GAPS_STATISTICS_H__
2 
+#define __GAPS_GAPS_STATISTICS_H__
3 
+
4 
+#include "GibbsSampler.h"
5 
+#include "Matrix.h"
6 
+
7 
+class GapsStatistics
8 
+{
9 
+private:
10 
+
11 
+    ColMatrix mAMeanMatrix;
12 
+    ColMatrix mAStdMatrix;
13 
+    RowMatrix mPMeanMatrix;
14 
+    RowMatrix mPStdMatrix;
15 
+
16 
+    unsigned mStatUpdates;
17 
+
18 
+public:
19 
+
20 
+    GapsStatistics(unsigned nRow, unsigned nCol, unsigned nFactor);
21 
+
22 
+    Rcpp::NumericMatrix AMean();
23 
+    Rcpp::NumericMatrix AStd();
24 
+    Rcpp::NumericMatrix PMean();
25 
+    Rcpp::NumericMatrix PStd();
26 
+
27 
+    void update(const AmplitudeGibbsSampler &ASampler,
28 
+        const PatternGibbsSampler &PSampler);
29 
+};
30 
+
31 
+#endif
0 32 
\ No newline at end of file
src/GibbsSampler.cpp
History View file @094c9c46c
... ... 
@@ -1,21 +1,20 @@
1 1 
 #include "GibbsSampler.h"
2 2
3 
-/*
4 3 
 AmplitudeGibbsSampler::AmplitudeGibbsSampler(const Rcpp::NumericMatrix &D,
5 
-const Rcpp::NumericMatrix &S, unsigned nFactor)
4 
+const Rcpp::NumericMatrix &S, unsigned nFactor, float alpha, float maxGibbsMass)
6 5 
     :
7 
-mMatrix(D.nrow(), nFactor), mDMatrix(D), mSMatrix(S), mAPMatrix(D.nrow(), D.ncol()),
8 
-mNumRows(D.nrow()), mNumCols(nFactor)
9 
-{}
6 
+mAMatrix(D.nrow(), nFactor), mDMatrix(D), mSMatrix(S),
7 
+mAPMatrix(D.nrow(), D.ncol()), mQueue(D.nrow() * nFactor, alpha),
8 
+mAnnealingTemp(0.f), mNumRows(D.nrow()), mNumCols(nFactor)
9 
+{
10 
+    mBinSize = std::numeric_limits<uint64_t>::max() / (mNumRows * mNumCols);
11 
+    uint64_t remain = std::numeric_limits<uint64_t>::max() % (mNumRows * mNumCols);
12 
+    mQueue.setDomainSize(std::numeric_limits<uint64_t>::max() - remain);
10 13
11 
-AmplitudeGibbsSampler::AmplitudeGibbsSampler(const Rcpp::NumericMatrix &D,
12 
-const Rcpp::NumericMatrix &S, unsigned nFactor, float alpha, float maxGibbsmass)
13 
-    :
14 
-mMatrix(D.nrow(), nFactor), mDMatrix(D), mSMatrix(S), mAPMatrix(D.nrow(), D.ncol()),
15 
-mQueue(alpha, D.ncol(), D.nrow() * nFactor), mMaxGibbsMass(maxGibbsmass),
16 
-mAnnealingTemp(0), mNumRows(D.nrow()), mNumCols(nFactor)
17 
-{}
18 
-*/
14 
+    float meanD = gaps::algo::mean(mDMatrix);
15 
+    mLambda = alpha * std::sqrt(nFactor / meanD);
16 
+    mMaxGibbsMass = maxGibbsMass / mLambda;
17 
+}
19 18
20 19 
 unsigned AmplitudeGibbsSampler::getRow(uint64_t pos) const
21 20 
 {
... ... 
@@ -29,39 +28,59 @@ unsigned AmplitudeGibbsSampler::getCol(uint64_t pos) const
29 28
30 29 
 bool AmplitudeGibbsSampler::canUseGibbs(unsigned row, unsigned col) const
31 30 
 {
32 
-    return !gaps::algo::isRowZero(*mOtherMatrix, col);
31 
+    return !gaps::algo::isRowZero(*mPMatrix, col);
33 32 
 }
34 33
35 34 
 bool AmplitudeGibbsSampler::canUseGibbs(unsigned r1, unsigned c1, unsigned r2, unsigned c2) const
36 35 
 {
37 
-    return !gaps::algo::isRowZero(*mOtherMatrix, c1)
38 
-        && !gaps::algo::isRowZero(*mOtherMatrix, c2);
36 
+    return !gaps::algo::isRowZero(*mPMatrix, c1)
37 
+        && !gaps::algo::isRowZero(*mPMatrix, c2);
39 38 
 }
40 39
41 40 
 void AmplitudeGibbsSampler::updateAPMatrix(unsigned row, unsigned col, float delta)
42 41 
 {
43 42 
     for (unsigned j = 0; j < mAPMatrix.nCol(); ++j)
44 43 
     {
45 
-        mAPMatrix(row,j) += delta * (*mOtherMatrix)(j,col);
44 
+        mAPMatrix(row,j) += delta * (*mPMatrix)(j,col);
46 45 
     }
47 46 
 }
48 47
49 
-/*
50 
-PatternGibbsSampler::PatternGibbsSampler(const Rcpp::NumericMatrix &D,
51 
-const Rcpp::NumericMatrix &S, unsigned nFactor)
52 
-    :
53 
-mMatrix(nFactor, D.ncol()), mDMatrix(D), mSMatrix(S), mAPMatrix(D.nrow(), D.ncol()),
54 
-mNumRows(nFactor), mNumCols(D.ncol())
55 
-{}
48 
+void AmplitudeGibbsSampler::sync(PatternGibbsSampler &sampler)
49 
+{
50 
+    mPMatrix = &(sampler.mPMatrix);
51 
+    mAPMatrix = sampler.mAPMatrix;
52 
+}
53 
+
54 
+float AmplitudeGibbsSampler::nAtoms() const
55 
+{
56 
+    return mDomain.size();
57 
+}
58 
+
59 
+void AmplitudeGibbsSampler::setAnnealingTemp(float temp)
60 
+{
61 
+    mAnnealingTemp = temp;
62 
+}
63 
+
64 
+float AmplitudeGibbsSampler::chi2() const
65 
+{
66 
+    return 2.f * gaps::algo::loglikelihood(mDMatrix, mSMatrix, mAPMatrix);
67 
+}
56 68
57 69 
 PatternGibbsSampler::PatternGibbsSampler(const Rcpp::NumericMatrix &D,
58 
-const Rcpp::NumericMatrix &S, unsigned nFactor, float alpha, float maxGibbsmass)
70 
+const Rcpp::NumericMatrix &S, unsigned nFactor, float alpha, float maxGibbsMass)
59 71 
     :
60 
-mMatrix(nFactor, D.ncol()), mDMatrix(D), mSMatrix(S), mAPMatrix(D.nrow(), D.ncol()),
61 
-mQueue(alpha, D.nrow(), D.ncol() * nFactor), mMaxGibbsMass(maxGibbsmass),
62 
-mAnnealingTemp(0), mNumRows(nFactor), mNumCols(D.ncol())
63 
-{}
64 
-*/
72 
+mPMatrix(nFactor, D.ncol()), mDMatrix(D), mSMatrix(S),
73 
+mAPMatrix(D.nrow(), D.ncol()), mQueue(nFactor * D.ncol(), alpha),
74 
+mAnnealingTemp(0.f), mNumRows(nFactor), mNumCols(D.nrow())
75 
+{
76 
+    mBinSize = std::numeric_limits<uint64_t>::max() / (mNumRows * mNumCols);
77 
+    uint64_t remain = std::numeric_limits<uint64_t>::max() % (mNumRows * mNumCols);
78 
+    mQueue.setDomainSize(std::numeric_limits<uint64_t>::max() - remain);
79 
+
80 
+    float meanD = gaps::algo::mean(mDMatrix);
81 
+    mLambda = alpha * std::sqrt(nFactor / meanD);
82 
+    mMaxGibbsMass = maxGibbsMass / mLambda;
83 
+}
65 84
66 85 
 unsigned PatternGibbsSampler::getRow(uint64_t pos) const
67 86 
 {
... ... 
@@ -75,19 +94,40 @@ unsigned PatternGibbsSampler::getCol(uint64_t pos) const
75 94
76 95 
 bool PatternGibbsSampler::canUseGibbs(unsigned row, unsigned col) const
77 96 
 {
78 
-    return !gaps::algo::isColZero(*mOtherMatrix, row);
97 
+    return !gaps::algo::isColZero(*mAMatrix, row);
79 98 
 }
80 99
81 100 
 bool PatternGibbsSampler::canUseGibbs(unsigned r1, unsigned c1, unsigned r2, unsigned c2) const
82 101 
 {
83 
-    return !gaps::algo::isColZero(*mOtherMatrix, r1)
84 
-        && !gaps::algo::isColZero(*mOtherMatrix, r2);
102 
+    return !gaps::algo::isColZero(*mAMatrix, r1)
103 
+        && !gaps::algo::isColZero(*mAMatrix, r2);
85 104 
 }
86 105
87 106 
 void PatternGibbsSampler::updateAPMatrix(unsigned row, unsigned col, float delta)
88 107 
 {
89 108 
     for (unsigned i = 0; i < mAPMatrix.nRow(); ++i)
90 109 
     {
91 
-        mAPMatrix(i,col) += delta * (*mOtherMatrix)(i,row);
110 
+        mAPMatrix(i,col) += delta * (*mAMatrix)(i,row);
92 111 
     }
93 112 
 }
113 
+
114 
+void PatternGibbsSampler::sync(AmplitudeGibbsSampler &sampler)
115 
+{
116 
+    mAMatrix = &(sampler.mAMatrix);
117 
+    mAPMatrix = sampler.mAPMatrix;
118 
+}
119 
+
120 
+float PatternGibbsSampler::nAtoms() const
121 
+{
122 
+    return mDomain.size();
123 
+}
124 
+
125 
+void PatternGibbsSampler::setAnnealingTemp(float temp)
126 
+{
127 
+    mAnnealingTemp = temp;
128 
+}
129 
+
130 
+float PatternGibbsSampler::chi2() const
131 
+{
132 
+    return 2.f * gaps::algo::loglikelihood(mDMatrix, mSMatrix, mAPMatrix);
133 
+}
94 134 
\ No newline at end of file
src/GibbsSampler.h
History View file @094c9c46c
... ... 
@@ -1,6 +1,8 @@
1 1 
 #ifndef __GAPS_GIBBS_SAMPLER_H__
2 2 
 #define __GAPS_GIBBS_SAMPLER_H__
3 3
4 
+#include "GapsAssert.h"
5 
+#include "GapsStatistics.h"
4 6 
 #include "Archive.h"
5 7 
 #include "Matrix.h"
6 8 
 #include "Random.h"
... ... 
@@ -8,26 +10,30 @@
8 10 
 #include "ProposalQueue.h"
9 11 
 #include "AtomicDomain.h"
10 12
11 
-// TODO have friend class for stats
12 
-// CRTP
13 
-template <class T, class MatA, class MatB>
14 
-class GibbsSampler
13 
+#include <Rcpp.h>
14 
+
15 
+// stats should be friend class
16 
+// CRTP not really needed here, free friend functions are clear enough,
17 
+// it would reduce some repeated code - but only a few lines and probably
18 
+// not worth the headache
19 
+// maybe convert to CRTP once everything is working cleanly
20 
+
21 
+class AmplitudeGibbsSampler;
22 
+class PatternGibbsSampler;
23 
+class GapsStatistics;
24 
+
25 
+class AmplitudeGibbsSampler
15 26 
 {
16 27 
 private:
17 28
18 
-    friend T; // prevent incorrect inheritance
19 
-    /*GibbsSampler(const Rcpp::NumericMatrix &D, const Rcpp::NumericMatrix &S);
20 
-    GibbsSampler(const Rcpp::NumericMatrix &D, const Rcpp::NumericMatrix &S,
21 
-        float alpha, float maxGibbsmass);*/
22 
-    GibbsSampler() {}
29 
+    friend PatternGibbsSampler;
30 
+    friend GapsStatistics;
23 31
24 
-protected:
25 
-
26 
-    MatA mMatrix;
27 
-    MatB* mOtherMatrix;
28 
-    MatB mDMatrix;
29 
-    MatB mSMatrix;
30 
-    MatB mAPMatrix;
32 
+    ColMatrix mAMatrix;
33 
+    RowMatrix *mPMatrix;
34 
+    RowMatrix mDMatrix;
35 
+    RowMatrix mSMatrix;
36 
+    RowMatrix mAPMatrix;
31 37
32 38 
     ProposalQueue mQueue;
33 39 
     AtomicDomain mDomain;
... ... 
@@ -35,61 +41,67 @@ protected:
35 41 
     float mLambda;
36 42 
     float mMaxGibbsMass;
37 43 
     float mAnnealingTemp;
38 
-
44 
+
39 45 
     unsigned mNumRows;
40 46 
     unsigned mNumCols;
41 47 
     unsigned mBinSize;
42 48
43 
-    T* impl();
44 
-
45 
-    void processProposal(const AtomicProposal &prop);
46 
-    void birth(uint64_t pos);
47 
-    void death(uint64_t pos, float mass);
48 
-    void move(uint64_t src, float mass, uint64_t dest);
49 
-    void exchange(uint64_t p1, float mass1, uint64_t p2, float mass2);
50 
-    float gibbsMass(unsigned row, unsigned col);
49 
+    unsigned getRow(uint64_t pos) const;
50 
+    unsigned getCol(uint64_t pos) const;
51 
+    bool canUseGibbs(unsigned row, unsigned col) const;
52 
+    bool canUseGibbs(unsigned r1, unsigned c1, unsigned r2, unsigned c2) const;
53 
+    void updateAPMatrix(unsigned row, unsigned col, float delta);
51 54
52 55 
 public:
53 56
54 
-    void update(unsigned nSteps);
55 
-    void syncAP(const MatA &otherAP);
56 
-    const MatB& APMatrix() const;
57 
-    void setAnnealingTemp(float temp);
57 
+    AmplitudeGibbsSampler(const Rcpp::NumericMatrix &D,
58 
+        const Rcpp::NumericMatrix &S, unsigned nFactor, float alpha=0.f,
59 
+        float maxGibbsMass=0.f);
58 60
59 
-    float chi2() const;
61 
+    void sync(PatternGibbsSampler &sampler);
60 62 
     float nAtoms() const;
63 
+    void setAnnealingTemp(float temp);
64 
+    float chi2() const;
61 65
62 
-    // serialization
63 
-    //friend Archive& operator<<(Archive &ar, GibbsSampler &sampler);
64 
-    //friend Archive& operator>>(Archive &ar, GibbsSampler &sampler);
65 
-};
66 
+    template <class Sampler>
67 
+    friend void update(Sampler&, unsigned);
66 68
67 
-class AmplitudeGibbsSampler : public GibbsSampler<AmplitudeGibbsSampler, ColMatrix, RowMatrix>
68 
-{
69 
-private:
69 
+    template <class Sampler>
70 
+    friend void processProposal(Sampler &sampler, const AtomicProposal &prop);
70 71
71 
-    friend GibbsSampler;
72 
+    template <class Sampler>
73 
+    friend void birth(Sampler &sampler, uint64_t pos);
72 74
73 
-    unsigned getRow(uint64_t pos) const;
74 
-    unsigned getCol(uint64_t pos) const;
75 
-    bool canUseGibbs(unsigned row, unsigned col) const;
76 
-    bool canUseGibbs(unsigned r1, unsigned c1, unsigned r2, unsigned c2) const;
77 
-    void updateAPMatrix(unsigned row, unsigned col, float delta);
78 
-
79 
-public:
75 
+    template <class Sampler>
76 
+    friend void death(Sampler &sampler, uint64_t pos, float mass);
80 77
81 
-    AmplitudeGibbsSampler(const Rcpp::NumericMatrix &D,
82 
-        const Rcpp::NumericMatrix &S, unsigned nFactor);
83 
-    AmplitudeGibbsSampler(const Rcpp::NumericMatrix &D,
84 
-        const Rcpp::NumericMatrix &S, unsigned nFactor, float alpha,
85 
-        float maxGibbsmass);
78 
+    template <class Sampler>
79 
+    friend void exchange(Sampler &sampler, uint64_t p1, float mass1, uint64_t p2, float mass2);
86 80 
 };
87 81
88 
-class PatternGibbsSampler : public GibbsSampler<PatternGibbsSampler, RowMatrix, ColMatrix>
82 
+class PatternGibbsSampler
89 83 
 {
90 84 
 private:
91 85
92 
-    friend GibbsSampler;
86 
+    friend AmplitudeGibbsSampler;
87 
+    friend GapsStatistics;
88 
+
89 
+    RowMatrix mPMatrix;
90 
+    ColMatrix *mAMatrix;
91 
+    ColMatrix mDMatrix;
92 
+    ColMatrix mSMatrix;
93 
+    ColMatrix mAPMatrix;
94 
+
95 
+    ProposalQueue mQueue;
96 
+    AtomicDomain mDomain;
97 
+
98 
+    float mLambda;
99 
+    float mMaxGibbsMass;
100 
+    float mAnnealingTemp;
101 
+
102 
+    unsigned mNumRows;
103 
+    unsigned mNumCols;
104 
+    unsigned mBinSize;
93 105
94 106 
     unsigned getRow(uint64_t pos) const;
95 107 
     unsigned getCol(uint64_t pos) const;
... ... 
@@ -100,28 +112,32 @@ private:
100 112 
 public:
101 113
102 114 
     PatternGibbsSampler(const Rcpp::NumericMatrix &D,
103 
-        const Rcpp::NumericMatrix &S, unsigned nFactor);
104 
-    PatternGibbsSampler(const Rcpp::NumericMatrix &D,
105 
-        const Rcpp::NumericMatrix &S, unsigned nFactor, float alpha,
106 
-        float maxGibbsmass);
107 
-};
115 
+        const Rcpp::NumericMatrix &S, unsigned nFactor, float alpha=0.f,
116 
+        float maxGibbsMass=0.f);
108 117
109 
-/*template <class T, class MatA, class MatB>
110 
-GibbsSampler<T,MatA,MatB>::GibbsSampler(const Rcpp::NumericMatrix &D,
111 
-const Rcpp::NumericMatrix &S, float alpha, float maxGibbsmass)
112 
-    :
113 
-mDMatrix(D), mSMatrix(S), mAPMatrix(D.nrow(), D.ncol()),
114 
-mMaxGibbsMass(maxGibbsmass), mAnnealingTemp(0.f)
115 
-{}*/
118 
+    void sync(AmplitudeGibbsSampler &sampler);
119 
+    float nAtoms() const;
120 
+    void setAnnealingTemp(float temp);
121 
+    float chi2() const;
116 122
117 
-template <class T, class MatA, class MatB>
118 
-T* GibbsSampler<T, MatA, MatB>::impl()
119 
-{
120 
-    return static_cast<T*>(this);
121 
-}
123 
+    template <class Sampler>
124 
+    friend void update(Sampler&, unsigned);
125 
+
126 
+    template <class Sampler>
127 
+    friend void processProposal(Sampler &sampler, const AtomicProposal &prop);
122 128
123 
-template <class T, class MatA, class MatB>
124 
-void GibbsSampler<T, MatA, MatB>::update(unsigned nSteps)
129 
+    template <class Sampler>
130 
+    friend void birth(Sampler &sampler, uint64_t pos);
131 
+
132 
+    template <class Sampler>
133 
+    friend void death(Sampler &sampler, uint64_t pos, float mass);
134 
+
135 
+    template <class Sampler>
136 
+    friend void exchange(Sampler &sampler, uint64_t p1, float mass1, uint64_t p2, float mass2);
137 
+};
138 
+
139 
+template <class Sampler>
140 
+void update(Sampler &sampler, unsigned nSteps)
125 141 
 {
126 142 
     unsigned n = 0;
127 143 
     while (n < nSteps)
... ... 
@@ -140,164 +156,208 @@ void GibbsSampler<T, MatA, MatB>::update(unsigned nSteps)
140 156 
         n += nJobs;
141 157 
         assert(n <= nSteps);
142 158 
         */
143 
-        mQueue.populate(mDomain, 1);
144 
-        processProposal(mQueue[0]);
145 
-        mQueue.clear();
159 
+        sampler.mQueue.populate(sampler.mDomain, 1);
160 
+        GAPS_ASSERT(sampler.mQueue.size() == 1);
161 
+        processProposal(sampler, sampler.mQueue[0]);
162 
+        sampler.mQueue.clear();
146 163 
         n++;
147 164 
     }
148 165 
 }
149 166
150 
-template <class T, class MatA, class MatB>
151 
-void GibbsSampler<T, MatA, MatB>::processProposal(const AtomicProposal &prop)
167 
+template <class Sampler>
168 
+void processProposal(Sampler &sampler, const AtomicProposal &prop)
152 169 
 {
170 
+    GAPS_ASSERT(prop.type == 'B' || prop.type == 'D' || prop.type == 'M' || prop.type == 'E');
153 171 
     switch (prop.type)
154 172 
     {
155 
-        case 'B': birth(prop.pos1); break;
156 
-        case 'D': death(prop.pos1, prop.mass1); break;
157 
-        case 'M': move(prop.pos1, prop.mass1, prop.pos2); break;
158 
-        case 'E': exchange(prop.pos1, prop.mass1, prop.pos2, prop.mass2); break;
173 
+        case 'B': birth(sampler, prop.pos1); break;
174 
+        case 'D': death(sampler, prop.pos1, prop.mass1); break;
175 
+        //case 'M': move(prop.pos1, prop.mass1, prop.pos2); break;
176 
+        case 'E': exchange(sampler, prop.pos1, prop.mass1, prop.pos2, prop.mass2); break;
159 177 
     }
160 178 
 }
161 179
162 
-template <class T, class MatA, class MatB>
163 
-void GibbsSampler<T, MatA, MatB>::birth(uint64_t pos)
180 
+template <class Sampler>
181 
+void birth(Sampler &sampler, uint64_t pos)
164 182 
 {
165 
-    unsigned row = impl()->getRow(pos);
166 
-    unsigned col = impl()->getCol(pos);
167 
-    //float mass = impl()->canUseGibbs(row, col) ? gibbsMass(row, col)
168 
-    //    : gaps::random::exponential(mLambda);
169 
-    float mass = gaps::random::exponential(mLambda);
170 
-
171 
-    mDomain.insert(pos, mass);
172 
-    mMatrix(row, col) += mass;
173 
-    impl()->updateAPMatrix(row, col, mass);
183 
+    unsigned row = sampler.getRow(pos);
184 
+    unsigned col = sampler.getCol(pos);
185 
+    float mass = gaps::random::exponential(sampler.mLambda);
186 
+    sampler.mDomain.insert(pos, mass);
174 187 
 }
175 188
176 
-template <class T, class MatA, class MatB>
177 
-void GibbsSampler<T, MatA, MatB>::death(uint64_t pos, float mass)
189 
+template <class Sampler>
190 
+void death(Sampler &sampler, uint64_t pos, float mass)
178 191 
 {
179 
-    unsigned row = impl()->getRow(pos);
180 
-    unsigned col = impl()->getCol(pos);
181 
-    mMatrix(row, col) += -mass;
182 
-    impl()->updateAPMatrix(row, col, -mass);
183 
-    mDomain.erase(pos);
184 
-    mQueue.acceptDeath();
185 
-
186 
-    /*float newMass = impl()->canUseGibbs(row, col) ? gibbsMass(row, col)
187 
-    : mass;
188 
-    float deltaLL = impl()->computeDeltaLL(row, col, newMass);
189 
-
190 
-    if (deltaLL * mAnnealingTemp >= std::log(gaps::random::uniform()))
191 
-    {
192 
-    float delta = mDomain.updateAtomMass(pos, newMass - mass);
193 
-    mMatrix(row, col) += delta;
194 
-    impl()->updateAPMatrix(row, col, delta);
195 
-    if (mDomain.at(pos))
196 
-    {
197 
-        mQueue.rejectDeath();
198 
-    }
199 
-    }
200 
-    else
201 
-    {
202 
-    mDomain.deleteAtom(pos);
203 
-    mQueue.maxAtoms--;
204 
-    }*/
192 
+    sampler.mQueue.rejectDeath();
205 193 
 }
206 194
207 
-template <class T, class MatA, class MatB>
208 
-void GibbsSampler<T, MatA, MatB>::move(uint64_t src, float mass, uint64_t dest)
195 
+template <class Sampler>
196 
+void exchange(Sampler &sampler, uint64_t p1, float mass1, uint64_t p2, float mass2)
209 197 
 {
210 
-    /*
211 
-    if (r1 == r2 && c1 == c2)
212 
-    {
213 
-    mDomain.deleteAtom(p1);
214 
-    mDomain.addAtom(p2, mass);
215 
-    }
216 
-    else
217 
-    {
218 
-    if (deltaLL * mAnnealingTemp >= std::log(gaps::random::uniform()))
219 
-    {
220 
-        mDomain.deleteAtom(p1);
221 
-        mDomain.addAtom(p2, mass);
222 
-        mMatrix(r1, c1) += -mass;
223 
-        mMatrix(r2, c2) += mass;
224 
-        impl()->updateAPMatrix(r1, c1, -mass);
225 
-        impl()->updateAPMatrix(r2, c2, mass);
226 
-    }
227 
-    }
228 
-    */
229 
-}
230 
-
231 
-template <class T, class MatA, class MatB>
232 
-void GibbsSampler<T, MatA, MatB>::exchange(uint64_t p1, float mass1, uint64_t p2, float mass2)
233 
-{
234 
-    /*
235 
-    float newMass = gaps::random::inverseGammaSample(0.f, mass1 + mass2, 2.f, 1.f / mLambda);
236 
-    float delta1 = mass1 > mass2 ? newMass - mass1 : mass2 - newMass;
237 
-    float delta2 = mass2 > mass1 ? newMass - mass2 : mass1 - newMass;
238 
-    unsigned r1 = impl()->getRow(p1);
239 
-    unsigned c1 = impl()->getCol(p1);
240 
-    unsigned r2 = impl()->getRow(p2);
241 
-    unsigned c2 = impl()->getCol(p2);
242 
-
243 
-    if (r1 == r2 && c1 == c2)
244 
-    {
245 
-    mDomain.updateAtomMass(p1, delta1);
246 
-    mDomain.updateAtomMass(p2, delta2);
247 
-    }
248 
-    else
249 
-    {
250 
-    // all the exchange code
251 
-    }
252 
-    */
253 
-    mQueue.rejectDeath();
254 
-}
255 
-
256 
-// don't return mass less than epsilon
257 
-template <class T, class MatA, class MatB>
258 
-float GibbsSampler<T, MatA, MatB>::gibbsMass(unsigned row, unsigned col)
259 
-{
260 
-/*
261 
-    AlphaParameters alpha = impl()->alphaParameters(row, col);
262 
-    alpha.s *= mAnnealingTemp / 2.f;
263 
-    alpha.su *= mAnnealingTemp / 2.f;
264 
-    float mean  = (2.f * alpha.su - mLambda) / (2.f * alpha.s);
265 
-    float sd = 1.f / std::sqrt(2.f * alpha.s);
266 
-
267 
-    float plower = gaps::random::p_norm(0.f, mean, sd);
268 
-
269 
-    float newMass = death ? mass : 0.f;
270 
-    if (plower < 1.f && alpha.s > 0.00001f)
271 
-    {
272 
-    newMass = gaps::random::inverseNormSample(plower, 1.f, mean, sd);
273 
-    }
274 
-    return std::max(0.f, std::min(newMax, mMaxGibbsMass));
275 
-*/
276 
-}
277 
-
278 
-template <class T, class MatA, class MatB>
279 
-void GibbsSampler<T, MatA, MatB>::syncAP(const MatA &otherAP)
280 
-{
281 
-    mAPMatrix = otherAP;
282 
-}
283 
-
284 
-template <class T, class MatA, class MatB>
285 
-const MatB& GibbsSampler<T, MatA, MatB>::APMatrix() const
286 
-{
287 
-    return mAPMatrix;
198 
+    sampler.mQueue.rejectDeath();
288 199 
 }
289 200
290 
-template <class T, class MatA, class MatB>
291 
-void GibbsSampler<T, MatA, MatB>::setAnnealingTemp(float temp)
292 
-{
293 
-    mAnnealingTemp = temp;
294 
-}
295 
-
296 
-template <class T, class MatA, class MatB>
297 
-float GibbsSampler<T, MatA, MatB>::chi2() const
298 
-{
299 
-    //return 2.f * gaps::algo::loglikelihood(mDMatrix, mSMatrix, mAPMatrix);
300 
-    return 0.f;
301 
-}
201 
+//  template <class T, class MatA, class MatB>
202 
+//  void GibbsSampler<T, MatA, MatB>::processProposal(const AtomicProposal &prop)
203 
+//  {
204 
+//      GAPS_ASSERT(prop.type == 'B' || prop.type == 'D' || prop.type == 'M' || prop.type == 'E');
205 
+//      switch (prop.type)
206 
+//      {
207 
+//          case 'B': birth(prop.pos1); break;
208 
+//          //case 'D': death(prop.pos1, prop.mass1); break;
209 
+//          //case 'M': move(prop.pos1, prop.mass1, prop.pos2); break;
210 
+//          //case 'E': exchange(prop.pos1, prop.mass1, prop.pos2, prop.mass2); break;
211 
+//      }
212 
+//  }
213 
+//
214 
+//  template <class T, class MatA, class MatB>
215 
+//  void GibbsSampler<T, MatA, MatB>::birth(uint64_t pos)
216 
+//  {
217 
+//      //GAPS_ASSERT(impl());
218 
+//      //unsigned row = impl()->getRow(pos);
219 
+//      //unsigned col = impl()->getCol(pos);
220 
+//      //float mass = impl()->canUseGibbs(row, col) ? gibbsMass(row, col)
221 
+//      //    : gaps::random::exponential(mLambda);
222 
+//      /*float mass = gaps::random::exponential(mLambda);
223 
+//
224 
+//      mDomain.insert(pos, mass);
225 
+//      mMatrix(row, col) += mass;
226 
+//      impl()->updateAPMatrix(row, col, mass);*/
227 
+//  }
228 
+//
229 
+//  template <class T, class MatA, class MatB>
230 
+//  void GibbsSampler<T, MatA, MatB>::death(uint64_t pos, float mass)
231 
+//  {
232 
+//      /*unsigned row = impl()->getRow(pos);
233 
+//      unsigned col = impl()->getCol(pos);
234 
+//      mMatrix(row, col) += -mass;
235 
+//      impl()->updateAPMatrix(row, col, -mass);
236 
+//      mDomain.erase(pos);
237 
+//      mQueue.acceptDeath();
238 
+//
239 
+//      float newMass = impl()->canUseGibbs(row, col) ? gibbsMass(row, col)
240 
+//      : mass;
241 
+//      float deltaLL = impl()->computeDeltaLL(row, col, newMass);
242 
+//
243 
+//      if (deltaLL * mAnnealingTemp >= std::log(gaps::random::uniform()))
244 
+//      {
245 
+//      float delta = mDomain.updateAtomMass(pos, newMass - mass);
246 
+//      mMatrix(row, col) += delta;
247 
+//      impl()->updateAPMatrix(row, col, delta);
248 
+//      if (mDomain.at(pos))
249 
+//      {
250 
+//          mQueue.rejectDeath();
251 
+//      }
252 
+//      }
253 
+//      else
254 
+//      {
255 
+//      mDomain.deleteAtom(pos);
256 
+//      mQueue.maxAtoms--;
257 
+//      }*/
258 
+//  }
259 
+//
260 
+//  template <class T, class MatA, class MatB>
261 
+//  void GibbsSampler<T, MatA, MatB>::move(uint64_t src, float mass, uint64_t dest)
262 
+//  {
263 
+//      /*
264 
+//      if (r1 == r2 && c1 == c2)
265 
+//      {
266 
+//      mDomain.deleteAtom(p1);
267 
+//      mDomain.addAtom(p2, mass);
268 
+//      }
269 
+//      else
270 
+//      {
271 
+//      if (deltaLL * mAnnealingTemp >= std::log(gaps::random::uniform()))
272 
+//      {
273 
+//          mDomain.deleteAtom(p1);
274 
+//          mDomain.addAtom(p2, mass);
275 
+//          mMatrix(r1, c1) += -mass;
276 
+//          mMatrix(r2, c2) += mass;
277 
+//          impl()->updateAPMatrix(r1, c1, -mass);
278 
+//          impl()->updateAPMatrix(r2, c2, mass);
279 
+//      }
280 
+//      }
281 
+//      */
282 
+//  }
283 
+//
284 
+//  template <class T, class MatA, class MatB>
285 
+//  void GibbsSampler<T, MatA, MatB>::exchange(uint64_t p1, float mass1, uint64_t p2, float mass2)
286 
+//  {
287 
+//      /*
288 
+//      float newMass = gaps::random::inverseGammaSample(0.f, mass1 + mass2, 2.f, 1.f / mLambda);
289 
+//      float delta1 = mass1 > mass2 ? newMass - mass1 : mass2 - newMass;
290 
+//      float delta2 = mass2 > mass1 ? newMass - mass2 : mass1 - newMass;
291 
+//      unsigned r1 = impl()->getRow(p1);
292 
+//      unsigned c1 = impl()->getCol(p1);
293 
+//      unsigned r2 = impl()->getRow(p2);
294 
+//      unsigned c2 = impl()->getCol(p2);
295 
+//
296 
+//      if (r1 == r2 && c1 == c2)
297 
+//      {
298 
+//      mDomain.updateAtomMass(p1, delta1);
299 
+//      mDomain.updateAtomMass(p2, delta2);
300 
+//      }
301 
+//      else
302 
+//      {
303 
+//      // all the exchange code
304 
+//      }
305 
+//      */
306 
+//      //mQueue.rejectDeath();
307 
+//  }
308 
+//
309 
+//  // don't return mass less than epsilon
310 
+//  template <class T, class MatA, class MatB>
311 
+//  float GibbsSampler<T, MatA, MatB>::gibbsMass(unsigned row, unsigned col)
312 
+//  {
313 
+//  /*
314 
+//      AlphaParameters alpha = impl()->alphaParameters(row, col);
315 
+//      alpha.s *= mAnnealingTemp / 2.f;
316 
+//      alpha.su *= mAnnealingTemp / 2.f;
317 
+//      float mean  = (2.f * alpha.su - mLambda) / (2.f * alpha.s);
318 
+//      float sd = 1.f / std::sqrt(2.f * alpha.s);
319 
+//
320 
+//      float plower = gaps::random::p_norm(0.f, mean, sd);
321 
+//
322 
+//      float newMass = death ? mass : 0.f;
323 
+//      if (plower < 1.f && alpha.s > 0.00001f)
324 
+//      {
325 
+//      newMass = gaps::random::inverseNormSample(plower, 1.f, mean, sd);
326 
+//      }
327 
+//      return std::max(0.f, std::min(newMax, mMaxGibbsMass));
328 
+//  */
329 
+//  }
330 
+//
331 
+//  template <class T, class MatA, class MatB>
332 
+//  void GibbsSampler<T, MatA, MatB>::syncAP(const MatA &otherAP)
333 
+//  {
334 
+//      mAPMatrix = otherAP;
335 
+//  }
336 
+//
337 
+//  template <class T, class MatA, class MatB>
338 
+//  const MatB& GibbsSampler<T, MatA, MatB>::APMatrix() const
339 
+//  {
340 
+//      return mAPMatrix;
341 
+//  }
342 
+//
343 
+//  template <class T, class MatA, class MatB>
344 
+//  void GibbsSampler<T, MatA, MatB>::setAnnealingTemp(float temp)
345 
+//  {
346 
+//      mAnnealingTemp = temp;
347 
+//  }
348 
+//
349 
+//  template <class Sampler>
350 
+//  float chi2(const Sampler &sampler)
351 
+//  {
352 
+//      //return 2.f * gaps::algo::loglikelihood(mDMatrix, mSMatrix, mAPMatrix);
353 
+//      return 0.f;
354 
+//  }
355 
+//
356 
+//  template <class T, class MatA, class MatB>
357 
+//  float GibbsSampler<T, MatA, MatB>::nAtoms() const
358 
+//  {
359 
+//      //return 2.f * gaps::algo::loglikelihood(mDMatrix, mSMatrix, mAPMatrix);
360 
+//      return 0.f;
361 
+//  }
302 362
303 363 
 #endif
304 364 
\ No newline at end of file
src/InternalState.h
History View file @094c9c46c
... ... 
@@ -16,6 +16,7 @@ enum GapsPhase
16 16 
     GAPS_SAMP
17 17 
 };
18 18
19 
+// maybe cleaner if this was a class with member functions? GapsRunner
19 20 
 struct GapsInternalState
20 21 
 {
21 22 
     Vector chi2VecEquil;
... ... 
@@ -58,7 +59,7 @@ struct GapsInternalState
58 59 
     GapsInternalState(const Rcpp::NumericMatrix &D,
59 60 
         const Rcpp::NumericMatrix &S, unsigned nF, unsigned nE, unsigned nEC,
60 61 
         unsigned nS, unsigned nOut, unsigned nSnap, float alphaA, float alphaP,
61 
-        float maxGibbmassA, float maxGibbmassP, int sd, bool msgs,
62 
+        float maxGibbsMassA, float maxGibbsMassP, int sd, bool msgs,
62 63 
         bool singleCellRNASeq, char whichMatrixFixed,
63 64 
         const Rcpp::NumericMatrix &FP, unsigned cptInterval)
64 65 
         //PumpThreshold pumpThreshold, unsigned numPumpSamples)
... ... 
@@ -69,8 +70,8 @@ struct GapsInternalState
69 70 
         nSnapshots(nSnap), nOutputs(nOut), messages(msgs), iter(0),
70 71 
         phase(GAPS_BURN), seed(sd), checkpointInterval(cptInterval),
71 72 
         nUpdatesA(0), nUpdatesP(0), //nPumpSamples(numPumpSamples),
72 
-        ASampler(D, S, nF, alphaA, maxGibbmassA),
73 
-        PSampler(D, S, nF, alphaP, maxGibbmassP)
73 
+        ASampler(D, S, nF, alphaA, maxGibbsMassA),
74 
+        PSampler(D, S, nF, alphaP, maxGibbsMassP)
74 75 
     {}
75 76
76 77 
     GapsInternalState(const Rcpp::NumericMatrix &D,
... ... 
@@ -84,6 +85,7 @@ struct GapsInternalState
84 85
85 86 
 inline Archive& operator<<(Archive &ar, GapsInternalState &state)
86 87 
 {
88 
+/*
87 89 
     ar << state.chi2VecEquil << state.nAtomsAEquil << state.nAtomsPEquil
88 90 
         << state.chi2VecSample << state.nAtomsASample << state.nAtomsPSample
89 91 
         << state.nIterA << state.nIterP << state.nEquil << state.nEquilCool
... ... 
@@ -91,11 +93,13 @@ inline Archive& operator<<(Archive &ar, GapsInternalState &state)
91 93 
         << state.iter << state.phase << state.seed << state.checkpointInterval
92 94 
         << state.nUpdatesA << state.nUpdatesP << state.nPumpSamples;
93 95 
         //<< state.sampler;
96 
+*/
94 97 
     return ar;
95 98 
 }
96 99
97 100 
 inline Archive& operator>>(Archive &ar, GapsInternalState &state)
98 101 
 {
102 
+/*
99 103 
     ar >> state.chi2VecEquil >> state.nAtomsAEquil >> state.nAtomsPEquil
100 104 
         >> state.chi2VecSample >> state.nAtomsASample >> state.nAtomsPSample
101 105 
         >> state.nIterA >> state.nIterP >> state.nEquil >> state.nEquilCool
... ... 
@@ -103,6 +107,7 @@ inline Archive& operator>>(Archive &ar, GapsInternalState &state)
103 107 
         >> state.iter >> state.phase >> state.seed >> state.checkpointInterval
104 108 
         >> state.nUpdatesA >> state.nUpdatesP >> state.nPumpSamples;
105 109 
         //>> state.sampler;
110 
+*/
106 111 
     return ar;
107 112 
 }
108 113
src/Makevars
History View file @094c9c46c
... ... 
@@ -1,4 +1,4 @@
1 
-PKG_CPPFLAGS = -DBOOST_MATH_PROMOTE_DOUBLE_POLICY=0 -DSIMD -msse4.1
1 
+PKG_CPPFLAGS = -DGAPS_DEBUG -DBOOST_MATH_PROMOTE_DOUBLE_POLICY=0 -DSIMD -msse4.1
2 2 
 OBJECTS =   Algorithms.o \
3 3 
             AtomicDomain.o \
4 4 
             Cogaps.o \
src/ProposalQueue.cpp
History View file @094c9c46c
... ... 
@@ -1,17 +1,36 @@
1 
+#include "GapsAssert.h"
1 2 
 #include "ProposalQueue.h"
2 3 
 #include "Random.h"
3 4
4 5 
 const uint64_t atomicEnd = std::numeric_limits<uint64_t>::max();
5 6 
 const double atomicSize = static_cast<double>(atomicEnd);
6 7
7 
-ProposalQueue::ProposalQueue(float alpha, uint64_t nIndices, uint64_t nBins)
8 
-    : mAlpha(alpha), mMinAtoms(0), mMaxAtoms(0), mNumIndices(nIndices), mNumBins(nBins)
9 
-{}
10 
-
11 8 
 void ProposalQueue::populate(const AtomicDomain &domain, unsigned limit)
12 9 
 {
13 
-    unsigned nIter = 0;
14 
-    while (nIter++ < limit && makeProposal(domain));
10 
+    //unsigned nIter = 0;
11 
+    //while (nIter++ < limit && makeProposal(domain));
12 
+    GAPS_ASSERT(makeProposal(domain));
13 
+    //Rprintf("%c\n", mQueue[0].type);
14 
+}
15 
+
16 
+void ProposalQueue::setNumBins(unsigned nBins)
17 
+{
18 
+    mNumBins = nBins;
19 
+}
20 
+
21 
+void ProposalQueue::setDomainSize(uint64_t size)
22 
+{
23 
+    mDomainSize = size;
24 
+}
25 
+
26 
+void ProposalQueue::setDimensionSize(unsigned size)
27 
+{
28 
+    mDimensionSize = size;
29 
+}
30 
+
31 
+void ProposalQueue::setAlpha(float alpha)
32 
+{
33 
+    mAlpha = alpha;
15 34 
 }
16 35
17 36 
 //void ProposalQueue::clear(unsigned n)
... ... 
@@ -29,9 +48,9 @@ unsigned ProposalQueue::size() const
29 48 
     return mQueue.size();
30 49 
 }
31 50
32 
-const AtomicProposal& ProposalQueue::operator[](unsigned n) const
51 
+const AtomicProposal& ProposalQueue::operator[](int n) const
33 52 
 {
34 
-    mQueue[mQueue.size() - 1 - n];
53 
+    return mQueue[mQueue.size() - 1 - n];
35 54 
 }
36 55
37 56 
 void ProposalQueue::acceptDeath()
... ... 
@@ -71,7 +90,7 @@ bool ProposalQueue::makeProposal(const AtomicDomain &domain)
71 90 
     float lowerBound = deleteProb(mMinAtoms);
72 91 
     float upperBound = deleteProb(mMaxAtoms);
73 92
74 
-    if (u1 < bdProb && u2 > upperBound)
93 
+    if (u1 < bdProb && u2 >= upperBound)
75 94 
     {
76 95 
         return birth(domain);
77 96 
     }
... ... 
@@ -86,7 +105,7 @@ bool ProposalQueue::makeProposal(const AtomicDomain &domain)
86 105 
     return false;
87 106 
 }
88 107
89 
-static bool isInVector(std::vector<unsigned> &vec, unsigned n)
108 
+static bool isInVector(const std::vector<unsigned> &vec, unsigned n)
90 109 
 {
91 110 
     return std::find(vec.begin(), vec.end(), n) != vec.end();
92 111 
 }
... ... 
@@ -94,13 +113,13 @@ static bool isInVector(std::vector<unsigned> &vec, unsigned n)
94 113 
 bool ProposalQueue::birth(const AtomicDomain &domain)
95 114 
 {
96 115 
     uint64_t pos = domain.randomFreePosition();
97 
-    if (isInVector(mUsedIndices, pos / mNumIndices))
98 
-    {
99 
-        return false; // matrix conflict - can't compute gibbs mass
100 
-    }
116 
+    //if (isInVector(mUsedIndices, pos / mDimensionSize))
117 
+    //{
118 
+    //    return false; // matrix conflict - can't compute gibbs mass
119 
+    //}
101 120 
     mQueue.push_back(AtomicProposal('B', pos));
102 
-    mUsedIndices.push_back(pos / mNumIndices);
103 
-    mUsedPositions.push_back(pos);
121 
+    //mUsedIndices.push_back(pos / mDimensionSize);
122 
+    //mUsedPositions.push_back(pos);
104 123 
     mMinAtoms++;
105 124 
     mMaxAtoms++;
106 125 
     return true;
... ... 
@@ -109,13 +128,13 @@ bool ProposalQueue::birth(const AtomicDomain &domain)
109 128 
 bool ProposalQueue::death(const AtomicDomain &domain)
110 129 
 {
111 130 
     Atom a = domain.randomAtom();
112 
-    if (isInVector(mUsedIndices, a.pos / mNumIndices))
113 
-    {
114 
-        return false; // matrix conflict - can't compute gibbs mass or deltaLL
115 
-    }
131 
+    //if (isInVector(mUsedIndices, a.pos / mDimensionSize))
132 
+    //{
133 
+    //    return false; // matrix conflict - can't compute gibbs mass or deltaLL
134 
+    //}
116 135 
     mQueue.push_back(AtomicProposal('D', a.pos, a.mass));
117 
-    mUsedIndices.push_back(a.pos / mNumIndices);
118 
-    mUsedPositions.push_back(a.pos);
136 
+    //mUsedIndices.push_back(a.pos / mDimensionSize);
137 
+    //mUsedPositions.push_back(a.pos);
119 138 
     mMinAtoms--;
120 139 
     return true;
121 140 
 }
... ... 
@@ -126,25 +145,25 @@ bool ProposalQueue::move(const AtomicDomain &domain)
126 145 
     uint64_t lbound = a.left ? a.left->pos : 0;
127 146 
     uint64_t rbound = a.right ? a.right->pos : atomicEnd - 1;
128 147
129 
-    for (unsigned i = 0; i < mUsedPositions.size(); ++i)
130 
-    {
131 
-        if (mUsedPositions[i] >= lbound && mUsedPositions[i] <= rbound)
132 
-        {
133 
-            return false; // atomic conflict - don't know neighbors
134 
-        }
135 
-    }
148 
+    //for (unsigned i = 0; i < mUsedPositions.size(); ++i)
149 
+    //{
150 
+    //    if (mUsedPositions[i] >= lbound && mUsedPositions[i] <= rbound)
151 
+    //    {
152 
+    //        return false; // atomic conflict - don't know neighbors
153 
+    //    }
154 
+    //}
136 155
137 156 
     uint64_t newLocation = gaps::random::uniform64(lbound, rbound);
138 
-    if (isInVector(mUsedIndices, a.pos / mNumIndices) || isInVector(mUsedIndices, newLocation / mNumIndices))
139 
-    {
140 
-        return false; // matrix conflict - can't compute deltaLL
141 
-    }
157 
+    //if (isInVector(mUsedIndices, a.pos / mDimensionSize) || isInVector(mUsedIndices, newLocation / mDimensionSize))
158 
+    //{
159 
+    //    return false; // matrix conflict - can't compute deltaLL
160 
+    //}
142 161
143 162 
     mQueue.push_back(AtomicProposal('M', a.pos, a.mass, newLocation));
144 
-    mUsedIndices.push_back(a.pos / mNumIndices);
145 
-    mUsedIndices.push_back(newLocation / mNumIndices);
146 
-    mUsedPositions.push_back(a.pos);
147 
-    mUsedPositions.push_back(newLocation);
163 
+    //mUsedIndices.push_back(a.pos / mDimensionSize);
164 
+    //mUsedIndices.push_back(newLocation / mDimensionSize);
165 
+    //mUsedPositions.push_back(a.pos);
166 
+    //mUsedPositions.push_back(newLocation);
148 167 
     return true;
149 168 
 }
150 169
... ... 
@@ -152,38 +171,38 @@ bool ProposalQueue::exchange(const AtomicDomain &domain)
152 171 
 {
153 172 
     Atom a1 = domain.randomAtom();
154 173
155 
-    if (a1.right) // has neighbor
156 
-    {
157 
-        for (unsigned i = 0; i < mUsedPositions.size(); ++i)
158 
-        {
159 
-            if (mUsedPositions[i] >= a1.pos && mUsedPositions[i] <= a1.right->pos)
160 
-            {
161 
-                return false; // atomic conflict - don't know right neighbor
162 
-            }
163 
-        }
164 
-    }
165 
-    else // exchange with first atom
166 
-    {
167 
-        for (unsigned i = 0; i < mUsedPositions.size(); ++i)
168 
-        {
169 
-            if (mUsedPositions[i] >= a1.pos || mUsedPositions[i] <= domain.front().pos)
170 
-            {
171 
-                return false; // atomic conflict - don't know right neighbor
172 
-            }
173 
-        }
174 
-    }
174 
+    //if (a1.right) // has neighbor
175 
+    //{
176 
+    //    for (unsigned i = 0; i < mUsedPositions.size(); ++i)
177 
+    //    {
178 
+    //        if (mUsedPositions[i] >= a1.pos && mUsedPositions[i] <= a1.right->pos)
179 
+    //        {
180 
+    //            return false; // atomic conflict - don't know right neighbor
181 
+    //        }
182 
+    //    }
183 
+    //}
184 
+    //else // exchange with first atom
185 
+    //{
186 
+    //    for (unsigned i = 0; i < mUsedPositions.size(); ++i)
187 
+    //    {
188 
+    //        if (mUsedPositions[i] >= a1.pos || mUsedPositions[i] <= domain.front().pos)
189 
+    //        {
190 
+    //            return false; // atomic conflict - don't know right neighbor
191 
+    //        }
192 
+    //    }
193 
+    //}
175 194
176 195 
     Atom a2 = a1.right ? *a1.right : domain.front();
177 
-    if (isInVector(mUsedIndices, a1.pos / mNumIndices) || isInVector(mUsedIndices, a2.pos / mNumIndices))
178 
-    {
179 
-        return false; // matrix conflict - can't compute gibbs mass or deltaLL
180 
-    }
196 
+    //if (isInVector(mUsedIndices, a1.pos / mDimensionSize) || isInVector(mUsedIndices, a2.pos / mDimensionSize))
197 
+    //{
198 
+    //    return false; // matrix conflict - can't compute gibbs mass or deltaLL
199 
+    //}
181 200
182 201 
     mQueue.push_back(AtomicProposal('E', a1.pos, a1.mass, a2.pos, a2.mass));
183 
-    mUsedIndices.push_back(a1.pos / mNumIndices);
184 
-    mUsedIndices.push_back(a2.pos / mNumIndices);
185 
-    mUsedPositions.push_back(a1.pos);
186 
-    mUsedPositions.push_back(a2.pos);
202 
+    //mUsedIndices.push_back(a1.pos / mDimensionSize);
203 
+    //mUsedIndices.push_back(a2.pos / mDimensionSize);
204 
+    //mUsedPositions.push_back(a1.pos);
205 
+    //mUsedPositions.push_back(a2.pos);
187 206 
     mMinAtoms--;
188 207 
     return true;
189 208 
 }
190 209 
\ No newline at end of file
src/ProposalQueue.h
History View file @094c9c46c
... ... 
@@ -31,13 +31,14 @@ private:
31 31 
     std::vector<AtomicProposal> mQueue; // not really a queue for now
32 32
33 33 
     std::vector<unsigned> mUsedIndices; // used rows/cols for A/P matrix
34 
-    std::vector<unsigned> mUsedPositions; // used positions is atomic domain
34 
+    std::vector<unsigned> mUsedPositions; // used positions in atomic domain
35 35
36 36 
     uint64_t mMinAtoms;
37 37 
     uint64_t mMaxAtoms;
38 38
39 39 
     uint64_t mNumBins;
40 
-    uint64_t mNumIndices;
40 
+    uint64_t mDimensionSize; // rows of A, cols of P
41 
+    uint64_t mDomainSize;
41 42
42 43 
     float mAlpha;
43 44
... ... 
@@ -50,14 +51,21 @@ private:
50 51
51 52 
 public:
52 53
53 
-    // constructor
54 
-    ProposalQueue(float alpha, uint64_t nIndices, uint64_t nBins);
54 
+    ProposalQueue(unsigned nBins, float alpha)
55 
+        : mMinAtoms(0), mMaxAtoms(0), mNumBins(nBins), mAlpha(alpha)
56 
+    {}
57 
+
58 
+    // set variables
59 
+    void setNumBins(unsigned nBins);
60 
+    void setDimensionSize(unsigned nIndices);
61 
+    void setDomainSize(uint64_t size);
62 
+    void setAlpha(float alpha);
55 63
56 64 
     // modify/access queue
57 65 
     void populate(const AtomicDomain &domain, unsigned limit);
58 66 
     void clear();
59 67 
     unsigned size() const;
60 
-    const AtomicProposal& operator[](unsigned n) const;
68 
+    const AtomicProposal& operator[](int n) const;
61 69
62 70 
     // update min/max atoms
63 71 
     void acceptDeath();