Developer documentation
Version 3.0.3-105-gd3941f44
random_threaded_loop.h
Go to the documentation of this file.
1/* Copyright (c) 2008-2022 the MRtrix3 contributors.
2 *
3 * This Source Code Form is subject to the terms of the Mozilla Public
4 * License, v. 2.0. If a copy of the MPL was not distributed with this
5 * file, You can obtain one at http://mozilla.org/MPL/2.0/.
6 *
7 * Covered Software is provided under this License on an "as is"
8 * basis, without warranty of any kind, either expressed, implied, or
9 * statutory, including, without limitation, warranties that the
10 * Covered Software is free of defects, merchantable, fit for a
11 * particular purpose or non-infringing.
12 * See the Mozilla Public License v. 2.0 for more details.
13 *
14 * For more details, see http://www.mrtrix.org/.
15 */
16
17#ifndef __algo_random_threaded_loop_h__
18#define __algo_random_threaded_loop_h__
19
20#include "debug.h"
21#include "exception.h"
22#include "algo/loop.h"
23#include "algo/iterator.h"
24#include "thread.h"
25#include "math/rng.h"
26#include <algorithm> // std::shuffle
27#include <random>
28// #include "algo/random_loop.h"
29
30namespace MR
31{
32
33 namespace {
34
35
36 template <int N, class Functor, class... ImageType>
37 struct RandomThreadedLoopRunInner
38 { MEMALIGN(RandomThreadedLoopRunInner<N,Functor,ImageType...>)
39 const vector<size_t>& outer_axes;
40 decltype (Loop (outer_axes)) loop;
41 typename std::remove_reference<Functor>::type func;
42 double density;
43 Math::RNG::Uniform<double> rng;
44 const vector<size_t> dims;
45 std::tuple<ImageType...> vox;
46
47 RandomThreadedLoopRunInner (const vector<size_t>& outer_axes, const vector<size_t>& inner_axes,
48 const Functor& functor, const double voxel_density, const vector<size_t> dimensions, ImageType&... voxels) :
49 outer_axes (outer_axes),
50 loop (Loop (inner_axes)),
51 func (functor),
52 density (voxel_density),
53 dims(dimensions),
54 vox (voxels...) { }
55
56 void operator() (const Iterator& pos) {
57 assign_pos_of (pos, outer_axes).to (vox);
58 for (auto i = unpack (loop, vox); i; ++i) {
59 // if (rng() >= density){
60 // DEBUG (str(pos) + " ...skipped inner");
61 // continue;
62 // }
63 // DEBUG (str(pos) + " ...used inner");
64 unpack (func, vox);
65 }
66 }
67 };
68
69
70 template <class Functor, class... ImageType>
71 struct RandomThreadedLoopRunInner<0,Functor,ImageType...>
72 { MEMALIGN(RandomThreadedLoopRunInner<0,Functor,ImageType...>)
73 const vector<size_t>& outer_axes;
74 const vector<size_t> inner;
75 decltype (Loop (outer_axes)) loop;
76 // Random_loop<Iterator, std::default_random_engine> random_loop;
77 typename std::remove_reference<Functor>::type func;
78 double density;
79 size_t max_cnt;
80 size_t cnt;
81 // Math::RNG::Uniform<double> rng;
82 std::default_random_engine random_engine;
83 vector<size_t> idx;
84 vector<size_t>::iterator it;
85 vector<size_t>::iterator stop;
86 const vector<size_t> dims;
87 // ImageType& image;
88 // RandomEngine& engine;
89 // size_t max_cnt;
90 // bool status;
91 // size_t cnt;
92
93
94 RandomThreadedLoopRunInner (const vector<size_t>& outer_axes, const vector<size_t>& inner_axes,
95 const Functor& functor, const double voxel_density, const vector<size_t>& dimensions, ImageType&... voxels) :
96 outer_axes (outer_axes),
97 inner (inner_axes),
98 loop (Loop (inner_axes)),
99 func (functor),
100 density (voxel_density),
101 dims (dimensions) {
102 assert (inner_axes.size() == 1);
103 // VAR(inner_axes);
104 // VAR(outer_axes);
105 Math::RNG rng;
106 typename std::default_random_engine::result_type seed = rng.get_seed();
107 random_engine = std::default_random_engine{static_cast<std::default_random_engine::result_type>(seed)};
108 idx = vector<size_t>(dims[inner_axes[0]]);
109 std::iota (std::begin(idx), std::end(idx), 0);
110 }
111
112 void operator() (Iterator& pos) {
113 std::shuffle (std::begin(idx), std::end(idx), random_engine);
114 cnt = 0;
115 it = std::begin(idx);
116 stop = std::end(idx);
117 max_cnt = //size_t (density * dims[inner[0]]);
118 (size_t) std::ceil (density * dims[inner[0]]);
119 // VAR(max_cnt);
120 for (size_t i = 0; i < max_cnt; ++i){
121 pos.index(inner[0]) = *it;
122 it++;
123 // for (auto i = loop (pos); i; ++i){
124 // if (rng() >= density){
125 // DEBUG (str(pos) + " ...skipped inner.");
126 // continue;
127 // }
128 // DEBUG (str(pos) + " ...used inner.");
129 // VAR(pos);
130 // VAR(cnt);
131 func (pos);
132 ++cnt;
133 }
134 pos.index(inner[0]) = dims[inner[0]];
135 }
136
137 };
138
139
140 template <class OuterLoopType>
141 struct RandomThreadedLoopRunOuter { MEMALIGN(RandomThreadedLoopRunOuter<OuterLoopType>)
142 Iterator iterator;
143 OuterLoopType outer_loop;
144 vector<size_t> inner_axes;
145
147 template <class Functor>
148 void run_outer (Functor&& functor, const double voxel_density, const vector<size_t>& dimensions)
149 {
150 if (Thread::threads_to_execute() == 0) {
151 for (auto i = outer_loop (iterator); i; ++i){
152 // std::cerr << "outer: " << str(iterator) << " " << voxel_density << " " << dimensions << std::endl;
153 functor (iterator);
154 }
155 return;
156 }
157
158 struct Shared { MEMALIGN(Shared)
159 Iterator& iterator;
160 decltype (outer_loop (iterator)) loop;
161 std::mutex mutex;
162 FORCE_INLINE bool next (Iterator& pos) {
163 std::lock_guard<std::mutex> lock (mutex);
164 if (loop) {
165 assign_pos_of (iterator, loop.axes).to (pos);
166 ++loop;
167 return true;
168 }
169 else return false;
170 }
171 } shared = { iterator, outer_loop (iterator) };
172
173 struct PerThread { MEMALIGN(PerThread)
174 Shared& shared;
175 typename std::remove_reference<Functor>::type func;
176 void execute () {
177 Iterator pos = shared.iterator;
178 while (shared.next (pos))
179 func (pos);
180 }
181 } loop_thread = { shared, functor };
182
183 Thread::run (Thread::multi (loop_thread), "loop threads").wait();
184 }
185
186
187
189 template <class Functor, class... ImageType>
190 void run (Functor&& functor, const double voxel_density, vector<size_t> dimensions, ImageType&&... vox)
191 {
192 RandomThreadedLoopRunInner<
193 sizeof...(ImageType),
194 typename std::remove_reference<Functor>::type,
195 typename std::remove_reference<ImageType>::type...
196 > loop_thread (outer_loop.axes, inner_axes, functor, voxel_density, dimensions, vox...);
197 run_outer (loop_thread, voxel_density, dimensions);
198 check_app_exit_code();
199 }
200
201 };
202 }
203
204
205
206
207
208
209
210 template <class HeaderType>
211 inline RandomThreadedLoopRunOuter<decltype(Loop(vector<size_t>()))> RandomThreadedLoop (
212 const HeaderType& source,
213 const vector<size_t>& outer_axes,
214 const vector<size_t>& inner_axes) {
215 return { source, Loop (outer_axes), inner_axes };
216 }
217
218
219 template <class HeaderType>
220 inline RandomThreadedLoopRunOuter<decltype(Loop(vector<size_t>()))> RandomThreadedLoop (
221 const HeaderType& source,
222 const vector<size_t>& axes,
223 size_t num_inner_axes = 1) {
224 return { source, Loop (get_outer_axes (axes, num_inner_axes)), get_inner_axes (axes, num_inner_axes) };
225 }
226
227 template <class HeaderType>
228 inline RandomThreadedLoopRunOuter<decltype(Loop(vector<size_t>()))> RandomThreadedLoop (
229 const HeaderType& source,
230 size_t from_axis = 0,
231 size_t to_axis = std::numeric_limits<size_t>::max(),
232 size_t num_inner_axes = 1) {
233 return { source,
234 Loop (get_outer_axes (source, num_inner_axes, from_axis, to_axis)),
235 get_inner_axes (source, num_inner_axes, from_axis, to_axis) };
236 }
237
238 template <class HeaderType>
239 inline RandomThreadedLoopRunOuter<decltype(Loop("", vector<size_t>()))> RandomThreadedLoop (
240 const std::string& progress_message,
241 const HeaderType& source,
242 const vector<size_t>& outer_axes,
243 const vector<size_t>& inner_axes) {
244 return { source, Loop (progress_message, outer_axes), inner_axes };
245 }
246
247 template <class HeaderType>
248 inline RandomThreadedLoopRunOuter<decltype(Loop("", vector<size_t>()))> RandomThreadedLoop (
249 const std::string& progress_message,
250 const HeaderType& source,
251 const vector<size_t>& axes,
252 size_t num_inner_axes = 1) {
253 return { source,
254 Loop (progress_message, get_outer_axes (axes, num_inner_axes)),
255 get_inner_axes (axes, num_inner_axes) };
256 }
257
258 template <class HeaderType>
259 inline RandomThreadedLoopRunOuter<decltype(Loop("", vector<size_t>()))> RandomThreadedLoop (
260 const std::string& progress_message,
261 const HeaderType& source,
262 size_t from_axis = 0,
263 size_t to_axis = std::numeric_limits<size_t>::max(),
264 size_t num_inner_axes = 1) {
265 return { source,
266 Loop (progress_message, get_outer_axes (source, num_inner_axes, from_axis, to_axis)),
267 get_inner_axes (source, num_inner_axes, from_axis, to_axis) };
268 }
269
270
273}
274
275#endif
276
277
278
run
void run()
MR::Math::RNG::get_seed
static std::mt19937::result_type get_seed()
Definition: rng.h:54
MR::Math::ceil
constexpr I ceil(const T x)
template function with cast to different type
Definition: math.h:86
MR::Loop
FORCE_INLINE LoopAlongAxes Loop()
Definition: loop.h:419
MR::Thread::run
__run< Functor >::type run(Functor &&functor, const std::string &name="unnamed")
Execute the functor's execute method in a separate thread.
Definition: thread.h:373
MR::Thread::threads_to_execute
size_t threads_to_execute()
MR::Thread::multi
__Multi< typename std::remove_reference< Functor >::type > multi(Functor &&functor, size_t nthreads=threads_to_execute())
used to request multiple threads of the corresponding functor
Definition: thread.h:285
MR::DWI::Tractography::rng
thread_local Math::RNG rng
thread-local, but globally accessible RNG to vastly simplify multi-threading
MR
Definition: base.h:24
MR::RandomThreadedLoop
RandomThreadedLoopRunOuter< decltype(Loop(vector< size_t >()))> RandomThreadedLoop(const HeaderType &source, const vector< size_t > &outer_axes, const vector< size_t > &inner_axes)
Definition: random_threaded_loop.h:211
MR::unpack
auto unpack(F &&f, T &&t) -> decltype(Unpack< ::std::tuple_size< typename ::std::decay< T >::type >::value >::unpack(::std::forward< F >(f), ::std::forward< T >(t)))
if t is a tuple of elements a..., invoke f(a...)
Definition: apply.h:91
MR::check_app_exit_code
void check_app_exit_code()
MEMALIGN
#define MEMALIGN(...)
Definition: types.h:185
FORCE_INLINE
#define FORCE_INLINE
Definition: types.h:156
functor
std::remove_reference< Functor >::type & functor
Definition: thread.h:215