2 Copyright (c) 2005-2019 Intel Corporation
4 Licensed under the Apache License, Version 2.0 (the "License");
5 you may not use this file except in compliance with the License.
6 You may obtain a copy of the License at
8 http://www.apache.org/licenses/LICENSE-2.0
10 Unless required by applicable law or agreed to in writing, software
11 distributed under the License is distributed on an "AS IS" BASIS,
12 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 See the License for the specific language governing permissions and
14 limitations under the License.
18 // Self-organizing map in TBB flow::graph
20 // This is an example of the use of cancellation in a graph. After a point in searching for
21 // the best match for an example, two examples are looked for simultaneously. When the
22 // earlier example is found and the update radius is determined, the affected searches
23 // for the subsequent example are cancelled, and after the update they are restarted.
24 // As the update radius shrinks fewer searches are cancelled, and by the last iterations
25 // virtually all the work done for the speculating example is useful.
27 // first, a simple implementation with only one example vector
30 // we will do a color map (the simple example.)
34 // for some number of iterations
35 // update radius r, weight of change L
37 // use graph to find BMU
38 // for each part of map within radius of BMU W
39 // update vector: W(t+1) = W(t) + w(dist)*L*(V - W(t))
44 #include "tbb/task_scheduler_init.h"
45 #include "tbb/flow_graph.h"
46 #include "tbb/blocked_range2d.h"
47 #include "tbb/tick_count.h"
48 #include "../../common/utility/utility.h"
54 static int xranges = 1;
55 static int yranges = 1;
56 static int xsize = -1;
57 static int ysize = -1;
59 static int global_i = 0;
60 static int speculation_start;
61 std::vector<int> function_node_execs;
62 static int xRangeMax = 3;
63 static int yRangeMax = 3;
64 static bool dont_speculate = false;
65 static search_result_type last_update;
67 class BMU_search_body {
69 subsquare_type my_square;
72 BMU_search_body(SOMap &_m, subsquare_type &_sq, int &fnt) : my_map(_m), my_square(_sq), fn_tally(fnt) { }
73 BMU_search_body( const BMU_search_body &other) : my_map(other.my_map), my_square(other.my_square), fn_tally(other.fn_tally) { }
74 search_result_type operator()(const SOM_element s) {
77 double min_dist = my_map.BMU_range(s, my_x, my_y, my_square);
78 ++fn_tally; // count how many times this function_node executed
79 return search_result_type(min_dist, my_x, my_y);
83 typedef function_node<SOM_element, search_result_type> search_node;
84 typedef broadcast_node<SOM_element> b_node;
85 typedef std::vector< search_node *> search_node_vector_type;
86 typedef std::vector< search_node_vector_type > search_node_array_type;
87 typedef std::vector< graph *> graph_vector_type;
88 typedef std::vector< graph_vector_type > graph_array_type;
90 #define SPECULATION_CNT 2
92 graph *g[SPECULATION_CNT]; // main graph; there should only be one per epoch
93 b_node *send_to[SPECULATION_CNT]; // broadcast node to send exemplar to all function_nodes
94 queue_node<search_result_type> *q[SPECULATION_CNT]; // queue for function nodes to put their results in
95 // each function_node should have its own graph
96 search_node_array_type* s_array[SPECULATION_CNT]; // 2d array of function nodes
97 graph_array_type* g_array[SPECULATION_CNT]; // 2d array of graphs
99 // build a set of SPECULATION_CNT graphs, each of which consists of a broadcast_node,
100 // xranges x yranges function_nodes, and one queue_node for output.
101 // once speculation starts, if i % SPECULATION_CNT is the current graph, (i+1) % SPECULATION_CNT
102 // is the first speculation, and so on.
104 build_BMU_graph(SOMap &map1) {
105 // build current graph
106 xsize = ((int)map1.size() + xranges - 1) / xranges;
107 ysize = ((int)map1[0].size() + yranges - 1) / yranges;
108 function_node_execs.clear();
109 function_node_execs.reserve(xranges*yranges+1);
110 for(int ii = 0; ii < xranges*yranges+1;++ii) function_node_execs.push_back(0);
112 for(int scnt = 0; scnt < SPECULATION_CNT; ++scnt) {
114 send_to[scnt] = new b_node(*(g[scnt])); // broadcast node to the function_nodes
115 q[scnt] = new queue_node<search_result_type>(*(g[scnt])); // output queue
117 // create the function_nodes, tie to the graph
118 s_array[scnt] = new search_node_array_type;
119 s_array[scnt]->reserve(xranges);
120 g_array[scnt] = new graph_array_type;
121 g_array[scnt]->reserve(xranges);
122 for(int i = 0; i < (int)map1.size(); i += xsize) {
123 int xindex = i / xsize;
124 s_array[scnt]->push_back(search_node_vector_type());
125 (*s_array[scnt])[xindex].reserve(yranges);
126 g_array[scnt]->push_back(graph_vector_type());
127 (*g_array[scnt])[xindex].reserve(yranges);
128 for( int j = 0; j < (int)map1[0].size(); j += ysize) {
129 int offset = (i/xsize)*yranges + (j / ysize);
130 int xmax = (i + xsize) > (int)map1.size() ? (int)map1.size() : i + xsize;
131 int ymax = (j + ysize) > (int)map1[0].size() ? (int)map1[0].size() : j + ysize;
132 subsquare_type sst(i,xmax,1,j,ymax,1);
133 BMU_search_body bb(map1,sst,function_node_execs[offset]);
134 graph *g_local = new graph;
135 search_node *s = new search_node(*g_local, serial, bb); // copies Body
136 (*g_array[scnt])[xindex].push_back(g_local);
137 (*s_array[scnt])[xindex].push_back(s);
138 make_edge(*(send_to[scnt]), *s); // broadcast_node -> function_node
139 make_edge(*s, *(q[scnt])); // function_node -> queue_node
145 // Wait for the 2D array of flow::graphs.
146 void wait_for_all_graphs(int cIndex) { // cIndex ranges over [0 .. SPECULATION_CNT - 1]
147 for(int x = 0; x < xranges; ++x) {
148 for(int y = 0; y < yranges; ++y) {
149 (*g_array[cIndex])[x][y]->wait_for_all();
155 destroy_BMU_graph() {
156 for(int scnt = 0; scnt < SPECULATION_CNT; ++scnt) {
157 for( int i = 0; i < (int)(*s_array[scnt]).size(); ++i ) {
158 for(int j = 0; j < (int)(*s_array[scnt])[i].size(); ++j) {
159 delete (*s_array[scnt])[i][j];
160 delete (*g_array[scnt])[i][j];
163 (*s_array[scnt]).clear();
164 delete s_array[scnt];
165 (*g_array[scnt]).clear();
166 delete g_array[scnt];
168 delete send_to[scnt];
173 void find_subrange_overlap(int const &xval, int const &yval, double const &radius, int &xlow, int &xhigh, int &ylow, int &yhigh) {
174 xlow = int((xval-radius)/xsize);
175 xhigh = int((xval+radius)/xsize);
176 ylow = int((yval-radius)/ysize);
177 yhigh = int((yval+radius)/ysize);
178 // circle may fall partly outside map
179 if(xlow < 0) xlow = 0;
180 if(xhigh >= xranges) xhigh = xranges - 1;
181 if(ylow < 0) ylow = 0;
182 if(yhigh >= yranges) yhigh = yranges - 1;
185 bool overlap( int &xval, int &yval, search_result_type &sr) {
186 int xlow, xhigh, ylow, yhigh;
187 find_subrange_overlap(get<XV>(sr), get<YV>(sr), get<RADIUS>(sr), xlow, xhigh, ylow, yhigh);
188 return xval >= xlow && xval <= xhigh && yval >= ylow && yval <= yhigh;
192 cancel_submaps(int &xval, int &yval, double &radius, int indx) {
197 find_subrange_overlap(xval, yval, radius, xlow, xhigh, ylow, yhigh);
198 for(int x = xlow; x <= xhigh; ++x) {
199 for(int y = ylow; y <= yhigh; ++y) {
200 (*g_array[indx])[x][y]->root_task()->cancel_group_execution();
206 restart_submaps(int &xval, int &yval, double &radius, int indx, SOM_element &vector) {
211 find_subrange_overlap(xval, yval, radius, xlow, xhigh, ylow, yhigh);
212 for(int x = xlow; x <= xhigh; ++x) {
213 for(int y = ylow; y <= yhigh; ++y) {
214 // have to reset the graph
215 (*g_array[indx])[x][y]->root_task()->context()->reset();
216 // and re-submit the exemplar for search.
217 (*s_array[indx])[x][y]->try_put(vector);
223 graph_BMU( int indx ) { // indx ranges over [0 .. SPECULATION_CNT -1]
224 wait_for_all_graphs(indx); // wait for the array of subgraphs
225 (g[indx])->wait_for_all();
226 std::vector<search_result_type> all_srs(xRangeMax*yRangeMax,search_result_type(DBL_MAX,-1,-1));
227 search_result_type sr;
228 search_result_type min_sr;
229 get<RADIUS>(min_sr) = DBL_MAX;
230 int result_count = 0;
231 while((q[indx])->try_get(sr)) {
233 // figure which submap this came from
234 int x = get<XV>(sr) / xsize;
235 int y = get<YV>(sr) / ysize;
236 int offset = x*yranges+y; // linearized subscript
237 all_srs[offset] = sr;
238 if(get<RADIUS>(sr) < get<RADIUS>(min_sr))
240 else if(get<RADIUS>(sr) == get<RADIUS>(min_sr)) {
241 if(get<XV>(sr) < get<XV>(min_sr)) {
244 else if((get<XV>(sr) == get<XV>(min_sr) &&
245 get<YV>(sr) < get<YV>(min_sr)))
255 void graph_teach(SOMap &map1, teaching_vector_type &in) {
256 build_BMU_graph(map1);
257 // normally the training would pick random exemplars to teach the SOM. We need
258 // the process to be reproducible, so we will pick the exemplars in order, [0, in.size())
260 for(int epoch = 0; epoch < nPasses; ++epoch) {
262 bool canceled_submaps = false;
263 int j = next_j; // try to make reproducible
264 next_j = (epoch+1) % in.size();
265 search_result_type min_sr;
266 if(epoch < speculation_start) {
267 (send_to[epoch%SPECULATION_CNT])->try_put(in[j]);
269 else if(epoch == speculation_start) {
270 (send_to[epoch%SPECULATION_CNT])->try_put(in[j]);
271 if(epoch < nPasses-1) {
272 (send_to[(epoch+1)%SPECULATION_CNT])->try_put(in[next_j]);
275 else if(epoch < nPasses - 1) {
276 (send_to[(epoch+1)%SPECULATION_CNT])->try_put(in[next_j]);
278 min_sr = graph_BMU(epoch % SPECULATION_CNT); //calls wait_for_all()
279 double min_distance = get<0>(min_sr);
280 double radius = max_radius * exp(-(double)epoch*radius_decay_rate);
281 double learning_rate = max_learning_rate * exp(-(double)epoch * learning_decay_rate);
282 if(epoch >= speculation_start && epoch < (nPasses - 1)) {
283 // have to cancel the affected submaps
284 cancel_submaps(get<XV>(min_sr), get<YV>(min_sr), radius, (epoch+1)%SPECULATION_CNT);
285 canceled_submaps = true;
287 map1.epoch_update(in[j], epoch, get<1>(min_sr), get<2>(min_sr), radius, learning_rate);
289 if(canceled_submaps) {
290 // do I have to wait for all the non-canceled speculative graph to complete first?
291 // yes, in case a canceled task was already executing.
292 wait_for_all_graphs((epoch+1) % SPECULATION_CNT); // wait for the array of subgraphs
293 restart_submaps(get<1>(min_sr), get<2>(min_sr), radius, (epoch+1)%SPECULATION_CNT, in[next_j]);
296 last_update = min_sr;
297 get<RADIUS>(last_update) = radius; // not smallest value, but range of effect
302 static const double serial_time_adjust = 1.25;
303 static double radius_fraction = 3.0;
306 main(int argc, char** argv) {
307 int l_speculation_start;
308 utility::thread_number_range threads(
309 task_scheduler_init::default_num_threads,
310 task_scheduler_init::default_num_threads() // run only the default number of threads if none specified
313 utility::parse_cli_arguments(argc,argv,
314 utility::cli_argument_pack()
315 //"-h" option for for displaying help is present implicitly
316 .positional_arg(threads,"n-of-threads","number of threads to use; a range of the form low[:high], where low and optional high are non-negative integers or 'auto' for the TBB default.")
317 // .positional_arg(InputFileName,"input-file","input file name")
318 // .positional_arg(OutputFileName,"output-file","output file name")
319 .positional_arg(radius_fraction, "radius-fraction","size of radius at which to start speculating")
320 .positional_arg(nPasses, "number-of-epochs","number of examples used in learning phase")
321 .arg(cancel_test, "cancel-test", "test for cancel signal while finding BMU")
322 .arg(extra_debug, "debug", "additional output")
323 .arg(dont_speculate,"nospeculate","don't speculate in SOM map teaching")
327 max_radius = (xMax < yMax) ? yMax / 2 : xMax / 2;
328 // need this value for the 1x1 timing below
329 radius_decay_rate = -(log(1.0/(double)max_radius) / (double)nPasses);
330 find_data_ranges(my_teaching, max_range, min_range );
332 printf( "Data range: ");
333 remark_SOM_element(min_range);
335 remark_SOM_element(max_range);
339 // find how much time is taken for the single function_node case.
340 // adjust nPasses so the 1x1 time is somewhere around serial_time_adjust seconds.
341 // make sure the example test runs for at least 0.5 second.
343 task_scheduler_init init(1);
344 SOMap map1(xMax,yMax);
345 speculation_start = nPasses + 1; // Don't speculate
349 map1.initialize(InitializeGradient, max_range, min_range);
350 tick_count t0 = tick_count::now();
351 graph_teach(map1, my_teaching);
352 tick_count t1 = tick_count::now();
353 double nSeconds = (t1-t0).seconds();
359 double size_adjust = sqrt(serial_time_adjust / nSeconds);
360 xMax = (int)((double)xMax * size_adjust);
361 yMax = (int)((double)yMax * size_adjust);
362 max_radius = (xMax < yMax) ? yMax / 2 : xMax / 2;
363 radius_decay_rate = log((double)max_radius) / (double)nPasses;
366 printf("original 1x1 case ran in %g seconds\n", nSeconds);
367 printf(" Size of table == %d x %d\n", xMax, yMax);
368 printf(" radius_decay_rate == %g\n", radius_decay_rate);
373 // the "max_radius" starts at 1/2*radius_fraction the table size. To start the speculation when the radius is
374 // 1 / n * the table size, the constant in the log below should be n / 2. so 2 == 1/4, 3 == 1/6th,
377 l_speculation_start = nPasses + 1;
378 if ( extra_debug )printf("speculation will not be done\n");
381 if(radius_fraction < 1.0 ) {
382 if ( extra_debug )printf("Warning: radius_fraction should be >= 1. Setting to 1.\n");
383 radius_fraction = 1.0;
385 l_speculation_start = (int)((double)nPasses * log(radius_fraction) / log((double)nPasses));
386 if ( extra_debug )printf( "We will start speculation at iteration %d\n", l_speculation_start );
388 double single_time; // for speedup calculations
389 for(int p = threads.first; p <= threads.last; ++p) {
390 task_scheduler_init init(p);
391 if ( extra_debug )printf( " -------------- Running with %d threads. ------------\n", p);
392 // run the SOM build for a series of subranges
393 for(xranges = 1; xranges <= xRangeMax; ++xranges) {
394 for(yranges = xranges; yranges <= yRangeMax; ++yranges) {
395 if(xranges == 1 && yranges == 1) {
396 // don't pointlessly speculate if we're only running one subrange.
397 speculation_start = nPasses + 1;
400 speculation_start = l_speculation_start;
402 SOMap map1(xMax, yMax);
403 map1.initialize(InitializeGradient, max_range, min_range);
405 if(extra_debug) printf( "Start learning for [%d,%d] ----------- \n", xranges,yranges);
406 tick_count t0 = tick_count::now();
407 graph_teach(map1, my_teaching);
408 tick_count t1 = tick_count::now();
410 if ( extra_debug )printf( "Done learning for [%d,%d], which took %g seconds ", xranges,yranges, (t1-t0).seconds());
411 if(xranges == 1 && yranges == 1) single_time = (t1-t0).seconds();
412 if ( extra_debug )printf( ": speedup == %g\n", single_time / (t1-t0).seconds());