xref: /oneTBB/examples/graph/som/som_graph.cpp (revision 6caecf96)
1 /*
2     Copyright (c) 2005-2021 Intel Corporation
3 
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
7 
8         http://www.apache.org/licenses/LICENSE-2.0
9 
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.
15 */
16 
17 //
18 // Self-organizing map in TBB flow::graph
19 //
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.
26 //
27 // first, a simple implementation with only one example vector
28 // at a time.
29 //
30 // we will do a color map (the simple example.)
31 //
32 //  graph algorithm
33 //
34 //       for some number of iterations
35 //           update radius r, weight of change L
36 //           for each example V
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))
40 
41 #define _MAIN_C_ 1
42 #include "som.hpp"
43 
44 #include "oneapi/tbb/flow_graph.h"
45 #include "oneapi/tbb/blocked_range2d.h"
46 #include "oneapi/tbb/tick_count.h"
47 #include "oneapi/tbb/task_arena.h"
48 
49 #include "common/utility/utility.hpp"
50 #include "common/utility/get_default_num_threads.hpp"
51 
52 #define RED   0
53 #define GREEN 1
54 #define BLUE  2
55 
56 static int xranges = 1;
57 static int yranges = 1;
58 static int xsize = -1;
59 static int ysize = -1;
60 
61 static int global_i = 0;
62 static int speculation_start;
63 
64 #if EXTRA_DEBUG
65 std::vector<int> cancel_count;
66 std::vector<int> extra_count;
67 std::vector<int> missing_count;
68 std::vector<int> canceled_before;
69 #endif
70 std::vector<int> function_node_execs;
71 static int xRangeMax = 3;
72 static int yRangeMax = 3;
73 static bool dont_speculate = false;
74 static search_result_type last_update;
75 
76 class BMU_search_body {
77     SOMap &my_map;
78     subsquare_type my_square;
79     int &fn_tally;
80 
81 public:
82     BMU_search_body(SOMap &_m, subsquare_type &_sq, int &fnt)
83             : my_map(_m),
84               my_square(_sq),
85               fn_tally(fnt) {}
86     BMU_search_body(const BMU_search_body &other)
87             : my_map(other.my_map),
88               my_square(other.my_square),
89               fn_tally(other.fn_tally) {}
90     search_result_type operator()(const SOM_element s) {
91         int my_x;
92         int my_y;
93         double min_dist = my_map.BMU_range(s, my_x, my_y, my_square);
94         ++fn_tally; // count how many times this function_node executed
95         return search_result_type(min_dist, my_x, my_y);
96     }
97 };
98 
99 typedef oneapi::tbb::flow::function_node<SOM_element, search_result_type> search_node;
100 typedef oneapi::tbb::flow::broadcast_node<SOM_element> b_node;
101 typedef std::vector<search_node *> search_node_vector_type;
102 typedef std::vector<search_node_vector_type> search_node_array_type;
103 typedef std::vector<oneapi::tbb::flow::graph *> graph_vector_type;
104 typedef std::vector<graph_vector_type> graph_array_type;
105 
106 #define SPECULATION_CNT 2
107 
108 oneapi::tbb::flow::graph *g[SPECULATION_CNT]; // main graph; there should only be one per epoch
109 b_node *send_to[SPECULATION_CNT]; // broadcast node to send exemplar to all function_nodes
110 oneapi::tbb::flow::queue_node<search_result_type>
111     *q[SPECULATION_CNT]; // queue for function nodes to put their results in
112 // each function_node should have its own graph
113 search_node_array_type *s_array[SPECULATION_CNT]; // 2d array of function nodes
114 graph_array_type *g_array[SPECULATION_CNT]; // 2d array of graphs
115 
116 // All graphs must locate in the same arena.
117 oneapi::tbb::flow::graph *construct_graph(oneapi::tbb::task_arena &ta) {
118     oneapi::tbb::flow::graph *result;
119     ta.execute([&result] {
120         result = new oneapi::tbb::flow::graph();
121     });
122     return result;
123 }
124 
125 // build a set of SPECULATION_CNT graphs, each of which consists of a broadcast_node,
126 //    xranges x yranges function_nodes, and one queue_node for output.
127 //    once speculation starts, if i % SPECULATION_CNT is the current graph, (i+1) % SPECULATION_CNT
128 //    is the first speculation, and so on.
129 void build_BMU_graph(SOMap &map1, oneapi::tbb::task_arena &ta) {
130     // build current graph
131     xsize = ((int)map1.size() + xranges - 1) / xranges;
132     ysize = ((int)map1[0].size() + yranges - 1) / yranges;
133     function_node_execs.clear();
134     function_node_execs.reserve(xranges * yranges + 1);
135     for (int i = 0; i < xranges * yranges + 1; ++i)
136         function_node_execs.push_back(0);
137 
138     for (int scnt = 0; scnt < SPECULATION_CNT; ++scnt) {
139         g[scnt] = construct_graph(ta);
140         send_to[scnt] = new b_node(*(g[scnt])); // broadcast node to the function_nodes
141         q[scnt] = new oneapi::tbb::flow::queue_node<search_result_type>(*(g[scnt])); // output queue
142 
143         // create the function_nodes, tie to the graph
144         s_array[scnt] = new search_node_array_type;
145         s_array[scnt]->reserve(xranges);
146         g_array[scnt] = new graph_array_type;
147         g_array[scnt]->reserve(xranges);
148         for (int i = 0; i < (int)map1.size(); i += xsize) {
149             int xindex = i / xsize;
150             s_array[scnt]->push_back(search_node_vector_type());
151 #if EXTRA_DEBUG
152             if (s_array[scnt]->size() != xindex + 1) {
153                 printf("Error; s_array[%d]->size() == %d, xindex== %d\n",
154                        scnt,
155                        (int)(s_array[scnt]->size()),
156                        xindex);
157             }
158 #endif
159             (*s_array[scnt])[xindex].reserve(yranges);
160             g_array[scnt]->push_back(graph_vector_type());
161             (*g_array[scnt])[xindex].reserve(yranges);
162             for (int j = 0; j < (int)map1[0].size(); j += ysize) {
163                 int offset = (i / xsize) * yranges + (j / ysize);
164                 int xmax = (i + xsize) > (int)map1.size() ? (int)map1.size() : i + xsize;
165                 int ymax = (j + ysize) > (int)map1[0].size() ? (int)map1[0].size() : j + ysize;
166                 subsquare_type sst(i, xmax, 1, j, ymax, 1);
167                 BMU_search_body bb(map1, sst, function_node_execs[offset]);
168                 oneapi::tbb::flow::graph *g_local = construct_graph(ta);
169                 search_node *s =
170                     new search_node(*g_local, oneapi::tbb::flow::serial, bb); // copies Body
171                 (*g_array[scnt])[xindex].push_back(g_local);
172                 (*s_array[scnt])[xindex].push_back(s);
173                 oneapi::tbb::flow::make_edge(*(send_to[scnt]),
174                                              *s); // broadcast_node -> function_node
175                 oneapi::tbb::flow::make_edge(*s, *(q[scnt])); // function_node -> queue_node
176             }
177         }
178     }
179 }
180 
181 // Wait for the 2D array of flow::graphs.
182 void wait_for_all_graphs(int cIndex) { // cIndex ranges over [0 .. SPECULATION_CNT - 1]
183     for (int x = 0; x < xranges; ++x) {
184         for (int y = 0; y < yranges; ++y) {
185             (*g_array[cIndex])[x][y]->wait_for_all();
186 #if EXTRA_DEBUG
187             __TBB_ASSERT(!(*g_array[cIndex])[x][y]->is_cancelled(),
188                          "wait_for_all() did not reset graph cancel");
189 #endif
190         }
191     }
192 }
193 
194 void destroy_BMU_graph() {
195     for (int scnt = 0; scnt < SPECULATION_CNT; ++scnt) {
196         for (int i = 0; i < (int)(*s_array[scnt]).size(); ++i) {
197             for (int j = 0; j < (int)(*s_array[scnt])[i].size(); ++j) {
198                 delete (*s_array[scnt])[i][j];
199                 delete (*g_array[scnt])[i][j];
200             }
201         }
202         (*s_array[scnt]).clear();
203         delete s_array[scnt];
204         (*g_array[scnt]).clear();
205         delete g_array[scnt];
206         delete q[scnt];
207         delete send_to[scnt];
208         delete g[scnt];
209     }
210 }
211 
212 void find_subrange_overlap(int const &xval,
213                            int const &yval,
214                            double const &radius,
215                            int &xlow,
216                            int &xhigh,
217                            int &ylow,
218                            int &yhigh) {
219     xlow = int((xval - radius) / xsize);
220     xhigh = int((xval + radius) / xsize);
221     ylow = int((yval - radius) / ysize);
222     yhigh = int((yval + radius) / ysize);
223     // circle may fall partly outside map
224     if (xlow < 0)
225         xlow = 0;
226     if (xhigh >= xranges)
227         xhigh = xranges - 1;
228     if (ylow < 0)
229         ylow = 0;
230     if (yhigh >= yranges)
231         yhigh = yranges - 1;
232 #if EXTRA_DEBUG
233     if (xlow >= xranges)
234         printf(" Error *** xlow == %d\n", xlow);
235     if (xhigh < 0)
236         printf("Error *** xhigh == %d\n", xhigh);
237     if (ylow >= yranges)
238         printf("Error *** ylow == %d\n", ylow);
239     if (yhigh < 0)
240         printf("Error *** yhigh == %d\n", yhigh);
241 #endif
242 }
243 
244 bool overlap(int &xval, int &yval, search_result_type &sr) {
245     int xlow, xhigh, ylow, yhigh;
246     find_subrange_overlap(
247         std::get<XV>(sr), std::get<YV>(sr), std::get<RADIUS>(sr), xlow, xhigh, ylow, yhigh);
248     return xval >= xlow && xval <= xhigh && yval >= ylow && yval <= yhigh;
249 }
250 
251 void cancel_submaps(int &xval, int &yval, double &radius, int indx) {
252     int xlow;
253     int xhigh;
254     int ylow;
255     int yhigh;
256     find_subrange_overlap(xval, yval, radius, xlow, xhigh, ylow, yhigh);
257     for (int x = xlow; x <= xhigh; ++x) {
258         for (int y = ylow; y <= yhigh; ++y) {
259             (*g_array[indx])[x][y]->cancel();
260         }
261     }
262 #if EXTRA_DEBUG
263     ++cancel_count[(xhigh - xlow + 1) * (yhigh - ylow + 1)];
264 #endif
265 }
266 
267 void restart_submaps(int &xval, int &yval, double &radius, int indx, SOM_element &vector) {
268     int xlow;
269     int xhigh;
270     int ylow;
271     int yhigh;
272     find_subrange_overlap(xval, yval, radius, xlow, xhigh, ylow, yhigh);
273     for (int x = xlow; x <= xhigh; ++x) {
274         for (int y = ylow; y <= yhigh; ++y) {
275             // have to reset the graph
276             (*g_array[indx])[x][y]->reset();
277             // and re-submit the exemplar for search.
278             (*s_array[indx])[x][y]->try_put(vector);
279         }
280     }
281 }
282 
283 search_result_type graph_BMU(int indx) { // indx ranges over [0 .. SPECULATION_CNT -1]
284     wait_for_all_graphs(indx); // wait for the array of subgraphs
285     (g[indx])->wait_for_all();
286     std::vector<search_result_type> all_srs(xRangeMax * yRangeMax,
287                                             search_result_type(DBL_MAX, -1, -1));
288 #if EXTRA_DEBUG
289     int extra_computations = 0;
290 #endif
291     search_result_type sr;
292     search_result_type min_sr;
293     std::get<RADIUS>(min_sr) = DBL_MAX;
294     int result_count = 0;
295     while ((q[indx])->try_get(sr)) {
296         ++result_count;
297         // figure which submap this came from
298         int x = std::get<XV>(sr) / xsize;
299         int y = std::get<YV>(sr) / ysize;
300 #if EXTRA_DEBUG
301         if (x < 0 || x >= xranges)
302             printf(" ###  x value out of range (%d)\n", x);
303         if (y < 0 || y >= yranges)
304             printf(" ###  y value out of range (%d)\n", y);
305 #endif
306         int offset = x * yranges + y; // linearized subscript
307 #if EXTRA_DEBUG
308         if (std::get<RADIUS>(all_srs[offset]) !=
309             DBL_MAX) { // we've already got a result from this subsquare
310             ++extra_computations;
311         }
312         else if (std::get<XV>(all_srs[offset]) != -1) {
313             if (extra_debug)
314                 printf("More than one cancellation of [%d,%d] iteration %d\n", x, y, global_i);
315         }
316 #endif
317         all_srs[offset] = sr;
318         if (std::get<RADIUS>(sr) < std::get<RADIUS>(min_sr))
319             min_sr = sr;
320         else if (std::get<RADIUS>(sr) == std::get<RADIUS>(min_sr)) {
321             if (std::get<XV>(sr) < std::get<XV>(min_sr)) {
322                 min_sr = sr;
323             }
324             else if ((std::get<XV>(sr) == std::get<XV>(min_sr) &&
325                       std::get<YV>(sr) < std::get<YV>(min_sr))) {
326                 min_sr = sr;
327             }
328         }
329     }
330 #if EXTRA_DEBUG
331     if (result_count != xranges * yranges + extra_computations) {
332         // we are missing at least one of the expected results.  Tally the missing values
333         for (int i = 0; i < xranges * yranges; ++i) {
334             if (std::get<RADIUS>(all_srs[i]) == DBL_MAX) {
335                 // i == x*yranges + y
336                 int xval = i / yranges;
337                 int yval = i % yranges;
338                 bool received_cancel_result = std::get<XV>(all_srs[i]) != -1;
339                 if (overlap(xval, yval, last_update)) {
340                     // we have previously canceled this subsquare.
341                     printf("No result for [%d,%d] which was canceled(%s)\n",
342                            xval,
343                            yval,
344                            received_cancel_result ? "T" : "F");
345                     ++canceled_before[i];
346                 }
347                 else {
348                     printf("No result for [%d,%d] which was not canceled(%s)\n",
349                            xval,
350                            yval,
351                            received_cancel_result ? "T" : "F");
352                 }
353                 ++missing_count[i];
354             }
355         }
356     }
357     if (extra_computations)
358         ++extra_count[extra_computations];
359 #endif
360     return min_sr;
361     // end of one epoch
362 }
363 
364 void graph_teach(SOMap &map1, teaching_vector_type &in, oneapi::tbb::task_arena &ta) {
365     build_BMU_graph(map1, ta);
366 #if EXTRA_DEBUG
367     cancel_count.clear();
368     extra_count.clear();
369     missing_count.clear();
370     canceled_before.clear();
371     cancel_count.reserve(xRangeMax * yRangeMax + 1);
372     extra_count.reserve(xRangeMax * yRangeMax + 1);
373     missing_count.reserve(xRangeMax * yRangeMax + 1);
374     canceled_before.reserve(xRangeMax * yRangeMax + 1);
375     for (int = 0; i < xRangeMax * yRangeMax + 1; ++i) {
376         cancel_count.push_back(0);
377         extra_count.push_back(0);
378         missing_count.push_back(0);
379         canceled_before.push_back(0);
380     }
381 #endif
382     // normally the training would pick random exemplars to teach the SOM.  We need
383     // the process to be reproducible, so we will pick the exemplars in order, [0, in.size())
384     int next_j = 0;
385     for (int epoch = 0; epoch < nPasses; ++epoch) {
386         global_i = epoch;
387         bool canceled_submaps = false;
388         int j = next_j; // try to make reproducible
389         next_j = (epoch + 1) % in.size();
390         search_result_type min_sr;
391         if (epoch < speculation_start) {
392             (send_to[epoch % SPECULATION_CNT])->try_put(in[j]);
393         }
394         else if (epoch == speculation_start) {
395             (send_to[epoch % SPECULATION_CNT])->try_put(in[j]);
396             if (epoch < nPasses - 1) {
397                 (send_to[(epoch + 1) % SPECULATION_CNT])->try_put(in[next_j]);
398             }
399         }
400         else if (epoch < nPasses - 1) {
401             (send_to[(epoch + 1) % SPECULATION_CNT])->try_put(in[next_j]);
402         }
403         min_sr = graph_BMU(epoch % SPECULATION_CNT); //calls wait_for_all()
404         double min_distance = std::get<0>(min_sr);
405         double radius = max_radius * exp(-(double)epoch * radius_decay_rate);
406         double learning_rate = max_learning_rate * exp(-(double)epoch * learning_decay_rate);
407         if (epoch >= speculation_start && epoch < (nPasses - 1)) {
408             // have to cancel the affected submaps
409             cancel_submaps(
410                 std::get<XV>(min_sr), std::get<YV>(min_sr), radius, (epoch + 1) % SPECULATION_CNT);
411             canceled_submaps = true;
412         }
413         map1.epoch_update(
414             in[j], epoch, std::get<1>(min_sr), std::get<2>(min_sr), radius, learning_rate);
415         ++global_i;
416         if (canceled_submaps) {
417             // do I have to wait for all the non-canceled speculative graph to complete first?
418             // yes, in case a canceled task was already executing.
419             wait_for_all_graphs((epoch + 1) % SPECULATION_CNT); // wait for the array of subgraphs
420             restart_submaps(std::get<1>(min_sr),
421                             std::get<2>(min_sr),
422                             radius,
423                             (epoch + 1) % SPECULATION_CNT,
424                             in[next_j]);
425         }
426 
427         last_update = min_sr;
428         std::get<RADIUS>(last_update) = radius; // not smallest value, but range of effect
429     }
430     destroy_BMU_graph();
431 }
432 
433 static const double serial_time_adjust = 1.25;
434 static double radius_fraction = 3.0;
435 
436 int main(int argc, char *argv[]) {
437     int l_speculation_start;
438     utility::thread_number_range threads(
439         utility::get_default_num_threads,
440         utility::
441             get_default_num_threads() // run only the default number of threads if none specified
442     );
443 
444     utility::parse_cli_arguments(
445         argc,
446         argv,
447         utility::cli_argument_pack()
448             //"-h" option for for displaying help is present implicitly
449             .positional_arg(
450                 threads,
451                 "n-of-threads",
452                 "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.")
453             // .positional_arg(InputFileName,"input-file","input file name")
454             // .positional_arg(OutputFileName,"output-file","output file name")
455             .positional_arg(
456                 radius_fraction, "radius-fraction", "size of radius at which to start speculating")
457             .positional_arg(
458                 nPasses, "number-of-epochs", "number of examples used in learning phase")
459             .arg(cancel_test, "cancel-test", "test for cancel signal while finding BMU")
460             .arg(extra_debug, "debug", "additional output")
461             .arg(dont_speculate, "nospeculate", "don't speculate in SOM map teaching"));
462 
463     readInputData();
464     max_radius = (xMax < yMax) ? yMax / 2 : xMax / 2;
465     // need this value for the 1x1 timing below
466     radius_decay_rate = -(log(1.0 / (double)max_radius) / (double)nPasses);
467     find_data_ranges(my_teaching, max_range, min_range);
468     if (extra_debug) {
469         printf("Data range: ");
470         remark_SOM_element(min_range);
471         printf(" to ");
472         remark_SOM_element(max_range);
473         printf("\n");
474     }
475 
476     // find how much time is taken for the single function_node case.
477     // adjust nPasses so the 1x1 time is somewhere around serial_time_adjust seconds.
478     // make sure the example test runs for at least 0.5 second.
479     for (;;) {
480         // Restrict max concurrency level via task_arena interface
481         oneapi::tbb::task_arena ta(1);
482         SOMap map1(xMax, yMax);
483         speculation_start = nPasses + 1; // Don't speculate
484 
485         xranges = 1;
486         yranges = 1;
487         map1.initialize(InitializeGradient, max_range, min_range);
488         oneapi::tbb::tick_count t0 = oneapi::tbb::tick_count::now();
489         graph_teach(map1, my_teaching, ta);
490         oneapi::tbb::tick_count t1 = oneapi::tbb::tick_count::now();
491         double nSeconds = (t1 - t0).seconds();
492         if (nSeconds < 0.5) {
493             xMax *= 2;
494             yMax *= 2;
495             continue;
496         }
497         double size_adjust = sqrt(serial_time_adjust / nSeconds);
498         xMax = (int)((double)xMax * size_adjust);
499         yMax = (int)((double)yMax * size_adjust);
500         max_radius = (xMax < yMax) ? yMax / 2 : xMax / 2;
501         radius_decay_rate = log((double)max_radius) / (double)nPasses;
502 
503         if (extra_debug) {
504             printf("original 1x1 case ran in %g seconds\n", nSeconds);
505             printf("   Size of table == %d x %d\n", xMax, yMax);
506             printf("   radius_decay_rate == %g\n", radius_decay_rate);
507         }
508         break;
509     }
510 
511     // the "max_radius" starts at 1/2*radius_fraction the table size.  To start the speculation when the radius is
512     // 1 / n * the table size, the constant in the log below should be n / 2.  so 2 == 1/4, 3 == 1/6th,
513     // et c.
514     if (dont_speculate) {
515         l_speculation_start = nPasses + 1;
516         if (extra_debug)
517             printf("speculation will not be done\n");
518     }
519     else {
520         if (radius_fraction < 1.0) {
521             if (extra_debug)
522                 printf("Warning: radius_fraction should be >= 1.  Setting to 1.\n");
523             radius_fraction = 1.0;
524         }
525         l_speculation_start = (int)((double)nPasses * log(radius_fraction) / log((double)nPasses));
526         if (extra_debug)
527             printf("We will start speculation at iteration %d\n", l_speculation_start);
528     }
529     double single_time; // for speedup calculations
530 #if EXTRA_DEBUG
531     // storage for the single-subrange answers, for comparing maps
532     std::vector<double> single_dist;
533     single_dist.reserve(my_teaching.size());
534     std::vector<int> single_xval;
535     single_xval.reserve(my_teaching.size());
536     std::vector<int> single_yval;
537     single_yval.reserve(my_teaching.size());
538 #endif
539     for (int p = threads.first; p <= threads.last; ++p) {
540         // Restrict max concurrency level via task_arena interface
541         oneapi::tbb::task_arena ta(p);
542         if (extra_debug)
543             printf(" -------------- Running with %d threads. ------------\n", p);
544         // run the SOM build for a series of subranges
545         for (xranges = 1; xranges <= xRangeMax; ++xranges) {
546             for (yranges = xranges; yranges <= yRangeMax; ++yranges) {
547                 if (xranges == 1 && yranges == 1) {
548                     // don't pointlessly speculate if we're only running one subrange.
549                     speculation_start = nPasses + 1;
550                 }
551                 else {
552                     speculation_start = l_speculation_start;
553                 }
554                 SOMap map1(xMax, yMax);
555                 map1.initialize(InitializeGradient, max_range, min_range);
556 
557                 if (extra_debug)
558                     printf("Start learning for [%d,%d] ----------- \n", xranges, yranges);
559                 oneapi::tbb::tick_count t0 = oneapi::tbb::tick_count::now();
560                 graph_teach(map1, my_teaching, ta);
561                 oneapi::tbb::tick_count t1 = oneapi::tbb::tick_count::now();
562 
563                 if (extra_debug)
564                     printf("Done learning for [%d,%d], which took %g seconds ",
565                            xranges,
566                            yranges,
567                            (t1 - t0).seconds());
568                 if (xranges == 1 && yranges == 1)
569                     single_time = (t1 - t0).seconds();
570                 if (extra_debug)
571                     printf(": speedup == %g\n", single_time / (t1 - t0).seconds());
572 
573 #if EXTRA_DEBUG
574                 if (extra_debug) {
575                     // number of times cancel was called, indexed by number of subranges canceled
576                     for (int i = 0; i < cancel_count.size(); ++i) {
577                         // only write output if we have a non-zero value.
578                         if (cancel_count[i] > 0) {
579                             int totalcnt = 0;
580                             printf("     cancellations: ");
581                             for (int j = 0; j < cancel_count.size(); ++j) {
582                                 if (cancel_count[j]) {
583                                     printf(" %d [%d]", j, cancel_count[j]);
584                                     totalcnt += cancel_count[j];
585                                 }
586                             }
587                             totalcnt += speculation_start;
588                             printf(" for a total of %d\n", totalcnt);
589                             break; // from for
590                         }
591                     }
592 
593                     // number of extra results (these occur when the subrange task starts before
594                     // cancel is received.)
595                     for (int i = 0; i < extra_count.size(); ++i) {
596                         if (extra_count[i] > 0) {
597                             int totalcnt = 0;
598                             printf("extra computations: ");
599                             for (int j = 0; j < extra_count.size(); ++j) {
600                                 if (extra_count[j]) {
601                                     printf(" %d[%d]", j, extra_count[j]);
602                                     totalcnt += extra_count[j];
603                                 }
604                             }
605                             totalcnt += speculation_start;
606                             printf(" for a total of %d\n", totalcnt);
607                             break; // from for
608                         }
609                     }
610 
611                     // here we count the number of times we looked for a particular subrange when fetching
612                     // the queue_node output and didn't find anything.  This may occur when a function_node
613                     // is "stuck" and doesn't process some number of exemplars.  function_node_execs is
614                     // a count of the number of times the corresponding function_node was executed (in
615                     // case the problem is dropped output in the queue_node.)
616                     for (int i = 0; i < missing_count.size(); ++i) {
617                         if (missing_count[i]) {
618                             int xval = i / yranges;
619                             int yval = i % yranges;
620                             printf(" f_node[%d,%d] missed %d values", xval, yval, missing_count[i]);
621                             if (canceled_before[i]) {
622                                 printf(" canceled_before == %d", canceled_before[i]);
623                             }
624                             printf(", fn_tally == %d\n", function_node_execs[i]);
625                         }
626                     }
627                 }
628 
629                 // check that output matches the 1x1 case
630                 for (int i = 0; i < my_teaching.size(); ++i) {
631                     int xdist;
632                     int ydist;
633                     double my_dist = map1.BMU(my_teaching[i], xdist, ydist);
634                     if (xranges == 1 && yranges == 1) {
635                         single_dist.push_back(my_dist);
636                         single_xval.push_back(xdist);
637                         single_yval.push_back(ydist);
638                     }
639                     else {
640                         if (single_dist[i] != my_dist || single_xval[i] != xdist ||
641                             single_yval[i] != ydist)
642                             printf(
643                                 "Error in output: expecting <%g, %d, %d>, but got <%g, %d, %d>\n",
644                                 single_dist[i],
645                                 single_xval[i],
646                                 single_yval[i],
647                                 my_dist,
648                                 xdist,
649                                 ydist);
650                     }
651                 }
652 #endif
653             } // yranges
654         } // xranges
655     } // #threads p
656     printf("done\n");
657     return 0;
658 }
659