Add index sampler
... | ... |
@@ -1,10 +1,84 @@ |
1 | 1 |
#include "GapsRunner.h" |
2 | 2 |
#include "math/SIMD.h" |
3 |
+#include "math/Random.h" |
|
4 |
+#include "GapsAssert.h" |
|
5 |
+#include <algorithm> |
|
3 | 6 |
|
4 | 7 |
#ifdef __GAPS_OPENMP__ |
5 | 8 |
#include <omp.h> |
6 | 9 |
#endif |
7 | 10 |
|
11 |
+// create "nSets" vectors where each vector contains a vector of indices in the |
|
12 |
+// range [0,n) |
|
13 |
+// see createGWCoGAPSSets.R - here we just sample indices, that function |
|
14 |
+// samples gene names as well |
|
15 |
+static std::vector< std::vector<unsigned> > sampleIndices(unsigned n, unsigned nSets) |
|
16 |
+{ |
|
17 |
+ unsigned setSize = (int) n / nSets; |
|
18 |
+ std::vector< std::vector<unsigned> > sampleIndices; |
|
19 |
+ std::vector<unsigned> toBeSampled; |
|
20 |
+ for (unsigned i = 1; i < n; ++i) |
|
21 |
+ { |
|
22 |
+ toBeSampled.push_back(i); |
|
23 |
+ } |
|
24 |
+ |
|
25 |
+ for (unsigned i = 0; i < (nSets - 1); ++i) |
|
26 |
+ { |
|
27 |
+ sampleIndices.push_back(gaps::random::sample(toBeSampled, setSize)); |
|
28 |
+ } |
|
29 |
+ |
|
30 |
+ GAPS_ASSERT(!toBeSampled.empty()); |
|
31 |
+ |
|
32 |
+ sampleIndices.push_back(toBeSampled); |
|
33 |
+ |
|
34 |
+ return sampleIndices; |
|
35 |
+ |
|
36 |
+ /* |
|
37 |
+ std::vector< std::vector<unsigned> > sampleIndices; |
|
38 |
+ std::vector<unsigned> sampled; |
|
39 |
+ |
|
40 |
+ for (unsigned i = 0; i < (n - 1) % nSets; ++i) |
|
41 |
+ { |
|
42 |
+ std::vector<unsigned> set; |
|
43 |
+ for (unsigned i = 0; i < ((unsigned) (n - 1) / nSets) + 1; ++i) |
|
44 |
+ { |
|
45 |
+ while(true) |
|
46 |
+ { |
|
47 |
+ unsigned sample = gaps::random::uniform64(1, n); |
|
48 |
+ if (find(sampled.begin(), sampled.end(), sample) == sampled.end()) |
|
49 |
+ { |
|
50 |
+ set.push_back(sample); |
|
51 |
+ sampled.push_back(sample); |
|
52 |
+ break; |
|
53 |
+ } |
|
54 |
+ } |
|
55 |
+ } |
|
56 |
+ sampleIndices.push_back(set); |
|
57 |
+ } |
|
58 |
+ |
|
59 |
+ for (unsigned i = (n - 1) % nSets; i < nSets; ++i) |
|
60 |
+ { |
|
61 |
+ std::vector<unsigned> set; |
|
62 |
+ for (unsigned i = 0; i < (unsigned) (n - 1) / nSets; ++i) |
|
63 |
+ { |
|
64 |
+ while(true) |
|
65 |
+ { |
|
66 |
+ unsigned sample = gaps::random::uniform64(1, n); |
|
67 |
+ if (find(sampled.begin(), sampled.end(), sample) == sampled.end()) |
|
68 |
+ { |
|
69 |
+ set.push_back(sample); |
|
70 |
+ sampled.push_back(sample); |
|
71 |
+ break; |
|
72 |
+ } |
|
73 |
+ } |
|
74 |
+ } |
|
75 |
+ sampleIndices.push_back(set); |
|
76 |
+ } |
|
77 |
+ |
|
78 |
+ return sampleIndices; |
|
79 |
+ */ |
|
80 |
+} |
|
81 |
+ |
|
8 | 82 |
GapsRunner::GapsRunner(const Rcpp::NumericMatrix &D, const Rcpp::NumericMatrix &S, |
9 | 83 |
unsigned nFactor, unsigned nEquil, unsigned nCool, unsigned nSample, |
10 | 84 |
unsigned nOutputs, unsigned nSnapshots, float alphaA, float alphaP, |
... | ... |
@@ -276,4 +350,4 @@ void GapsRunner::makeCheckpointIfNeeded() |
276 | 350 |
createCheckpoint(); |
277 | 351 |
mLastCheckpoint = bpt_now(); |
278 | 352 |
} |
279 |
-} |
|
280 | 353 |
\ No newline at end of file |
354 |
+} |
... | ... |
@@ -21,6 +21,7 @@ |
21 | 21 |
#include <boost/random/uniform_real_distribution.hpp> |
22 | 22 |
|
23 | 23 |
#include <stdint.h> |
24 |
+#include <algorithm> |
|
24 | 25 |
|
25 | 26 |
#ifdef __GAPS_OPENMP__ |
26 | 27 |
#include <omp.h> |
... | ... |
@@ -178,4 +179,37 @@ float gaps::random::inverseGammaSample(float a, float b, float mean, float sd) |
178 | 179 |
u = gaps::random::uniform(a, b); |
179 | 180 |
} |
180 | 181 |
return gaps::random::q_gamma(u, mean, sd); |
181 |
-} |
|
182 | 182 |
\ No newline at end of file |
183 |
+} |
|
184 |
+ |
|
185 |
+std::vector<unsigned> sample(std::vector<unsigned> &elements, unsigned n) { |
|
186 |
+ std::vector<unsigned> sampleVect; |
|
187 |
+ for (unsigned i = 0; i < n; ++i) |
|
188 |
+ { |
|
189 |
+ GAPS_ASSERT(n <= elements.size()); |
|
190 |
+ unsigned sampleIndex = gaps::random::uniform64(0, elements.size()); |
|
191 |
+ sampleVect.push_back(elements.at(sampleIndex)); |
|
192 |
+ |
|
193 |
+ elements[sampleIndex] = elements.back(); |
|
194 |
+ elements.pop_back(); |
|
195 |
+ } |
|
196 |
+ return sampleVect; |
|
197 |
+ |
|
198 |
+ /* |
|
199 |
+ std::vector<unsigned> sampleVect; |
|
200 |
+ std::vector<unsigned> sampledIndices; |
|
201 |
+ for (unsigned i = 0; i < n; ++i) |
|
202 |
+ { |
|
203 |
+ while(true) |
|
204 |
+ { |
|
205 |
+ unsigned sampleIndex = gaps::random::uniform64(0, elements.size()); |
|
206 |
+ if (find(sampledIndices.begin(), sampledIndices.end(), sampleIndex) == sampledIndices.end()) |
|
207 |
+ { |
|
208 |
+ sampleVect.push_back(elements.at(sampleIndex)); |
|
209 |
+ sampledIndices.push_back(sampleIndex); |
|
210 |
+ break; |
|
211 |
+ } |
|
212 |
+ } |
|
213 |
+ } |
|
214 |
+ return sampleVect; |
|
215 |
+ */ |
|
216 |
+} |
... | ... |
@@ -34,7 +34,9 @@ namespace gaps |
34 | 34 |
|
35 | 35 |
void save(Archive &ar); |
36 | 36 |
void load(Archive &ar); |
37 |
+ |
|
38 |
+ std::vector<unsigned> sample(const std::vector<unsigned> &elements, unsigned n); |
|
37 | 39 |
} // namespace random |
38 | 40 |
} // namespace gaps |
39 | 41 |
|
40 |
-#endif |
|
41 | 42 |
\ No newline at end of file |
43 |
+#endif |