1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
|
/* Copyright (c) 2008-2022 the MRtrix3 contributors.
*
* This Source Code Form is subject to the terms of the Mozilla Public
* License, v. 2.0. If a copy of the MPL was not distributed with this
* file, You can obtain one at http://mozilla.org/MPL/2.0/.
*
* Covered Software is provided under this License on an "as is"
* basis, without warranty of any kind, either expressed, implied, or
* statutory, including, without limitation, warranties that the
* Covered Software is free of defects, merchantable, fit for a
* particular purpose or non-infringing.
* See the Mozilla Public License v. 2.0 for more details.
*
* For more details, see http://www.mrtrix.org/.
*/
#ifndef __math_stats_shuffle_h__
#define __math_stats_shuffle_h__
#include "app.h"
#include "progressbar.h"
#include "types.h"
#include "misc/bitset.h"
#include "math/stats/typedefs.h"
#define DEFAULT_NUMBER_SHUFFLES 5000
#define DEFAULT_NUMBER_SHUFFLES_NONSTATIONARITY 5000
namespace MR
{
namespace Math
{
namespace Stats
{
// Generic command-line options:
// - Set nature of errors
// - Set number of shuffles (actual & nonstationarity correction)
// - Import permutations (actual & nonstationarity correction)
// - (future) Set exchangeability blocks
extern const char* error_types[];
App::OptionGroup shuffle_options (const bool include_nonstationarity, const default_type default_skew = 1.0);
class Shuffle
{ NOMEMALIGN
public:
size_t index;
matrix_type data;
};
class Shuffler
{ NOMEMALIGN
public:
typedef vector<size_t> PermuteLabels;
enum class error_t { EE, ISE, BOTH };
// First version reads command-line options in order to determine parameters prior to running initialise();
// second and third versions more-or-less call initialise() directly
Shuffler (const size_t num_rows,
const bool is_nonstationarity,
const std::string msg = "");
Shuffler (const size_t num_rows,
const size_t num_shuffles,
const error_t error_types,
const bool is_nonstationarity,
const std::string msg = "");
Shuffler (const size_t num_rows,
const size_t num_shuffles,
const error_t error_types,
const bool is_nonstationarity,
const index_array_type& eb_within,
const index_array_type& eb_whole,
const std::string msg = "");
// Don't store the full set of shuffling matrices;
// generate each as it is required, based on the more compressed representations
bool operator() (Shuffle& output);
size_t size() const { return nshuffles; }
// Go back to the first permutation
void reset();
private:
const size_t rows;
vector<PermuteLabels> permutations;
vector<BitSet> signflips;
size_t nshuffles, counter;
std::unique_ptr<ProgressBar> progress;
void initialise (const error_t error_types,
const bool nshuffles_explicit,
const bool is_nonstationarity,
const index_array_type& eb_within,
const index_array_type& eb_whole);
// For exchangeability blocks (either within or whole)
index_array_type load_blocks (const std::string& filename, const bool equal_sizes);
// For generating unique permutations
bool is_duplicate (const PermuteLabels&, const PermuteLabels&) const;
bool is_duplicate (const PermuteLabels&) const;
// Note that this function does not take into account identical rows and therefore generated
// permutations are not guaranteed to be unique wrt the computed test statistic.
// Providing the number of rows is large then the likelihood of generating duplicates is low.
void generate_random_permutations (const size_t num_perms,
const size_t num_rows,
const index_array_type& eb_within,
const index_array_type& eb_whole,
const bool include_default,
const bool permit_duplicates);
void generate_all_permutations (const size_t num_rows,
const index_array_type& eb_within,
const index_array_type& eb_whole);
void load_permutations (const std::string& filename);
// Similar functions required for sign-flipping
bool is_duplicate (const BitSet&) const;
void generate_random_signflips (const size_t num_signflips,
const size_t num_rows,
const index_array_type& blocks,
const bool include_default,
const bool permit_duplicates);
void generate_all_signflips (const size_t num_rows,
const index_array_type& blocks);
vector<vector<size_t>> indices2blocks (const index_array_type&) const;
};
}
}
}
#endif
|
