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 // we will do a color map (the simple example.)
21 //
22 // serial algorithm
23 //
24 // initialize map with vectors (could be random, gradient, or something else)
25 // for some number of iterations
26 // update radius r, weight of change L
27 // for each example V
28 // find the best matching unit
29 // for each part of map within radius of BMU W
30 // update vector: W(t+1) = W(t) + w(dist)*L*(V - W(t))
31
32 #include "oneapi/tbb/task_group.h"
33
34 #include "som.hpp"
35
operator <<(std::ostream & out,const SOM_element & s)36 std::ostream &operator<<(std::ostream &out, const SOM_element &s) {
37 out << "(";
38 for (int i = 0; i < (int)s.w.size(); ++i) {
39 out << s.w[i];
40 if (i < (int)s.w.size() - 1) {
41 out << ",";
42 }
43 }
44 out << ")";
45 return out;
46 }
47
remark_SOM_element(const SOM_element & s)48 void remark_SOM_element(const SOM_element &s) {
49 printf("(");
50 for (int i = 0; i < (int)s.w.size(); ++i) {
51 printf("%g", s.w[i]);
52 if (i < (int)s.w.size() - 1) {
53 printf(",");
54 }
55 }
56 printf(")");
57 }
58
operator <<(std::ostream & out,const search_result_type & s)59 std::ostream &operator<<(std::ostream &out, const search_result_type &s) {
60 out << "<";
61 out << std::get<RADIUS>(s);
62 out << ", " << std::get<XV>(s);
63 out << ", ";
64 out << std::get<YV>(s);
65 out << ">";
66 return out;
67 }
68
remark_search_result_type(const search_result_type & s)69 void remark_search_result_type(const search_result_type &s) {
70 printf("<%g,%d,%d>", std::get<RADIUS>(s), std::get<XV>(s), std::get<YV>(s));
71 }
72
randval(double lowlimit,double highlimit)73 double randval(double lowlimit, double highlimit) {
74 return double(rand()) / double(RAND_MAX) * (highlimit - lowlimit) + lowlimit;
75 }
76
find_data_ranges(teaching_vector_type & teaching,SOM_element & max_range,SOM_element & min_range)77 void find_data_ranges(teaching_vector_type &teaching,
78 SOM_element &max_range,
79 SOM_element &min_range) {
80 if (teaching.size() == 0)
81 return;
82 max_range = min_range = teaching[0];
83 for (int i = 1; i < (int)teaching.size(); ++i) {
84 max_range.elementwise_max(teaching[i]);
85 min_range.elementwise_min(teaching[i]);
86 }
87 }
88
add_fraction_of_difference(SOM_element & to,SOM_element const & from,double frac)89 void add_fraction_of_difference(SOM_element &to, SOM_element const &from, double frac) {
90 for (int i = 0; i < (int)from.size(); ++i) {
91 to[i] += frac * (from[i] - to[i]);
92 }
93 }
94
distance_squared(SOM_element x,SOM_element y)95 double distance_squared(SOM_element x, SOM_element y) {
96 double rval = 0.0;
97 for (int i = 0; i < (int)x.size(); ++i) {
98 double diff = x[i] - y[i];
99 rval += diff * diff;
100 }
101 return rval;
102 }
103
initialize(InitializeType it,SOM_element & max_range,SOM_element & min_range)104 void SOMap::initialize(InitializeType it, SOM_element &max_range, SOM_element &min_range) {
105 for (int x = 0; x < xMax; ++x) {
106 for (int y = 0; y < yMax; ++y) {
107 for (int i = 0; i < (int)max_range.size(); ++i) {
108 if (it == InitializeRandom) {
109 my_map[x][y][i] = (randval(min_range[i], max_range[i]));
110 }
111 else if (it == InitializeGradient) {
112 my_map[x][y][i] =
113 ((double)(x + y) / (xMax + yMax) * (max_range[i] - min_range[i]) +
114 min_range[i]);
115 }
116 }
117 }
118 }
119 }
120
121 // subsquare [low,high)
BMU_range(const SOM_element & s,int & xval,int & yval,subsquare_type & r)122 double SOMap::BMU_range(const SOM_element &s, int &xval, int &yval, subsquare_type &r) {
123 double min_distance_squared = DBL_MAX;
124 int min_x = -1;
125 int min_y = -1;
126 for (int x = r.rows().begin(); x != r.rows().end(); ++x) {
127 for (int y = r.cols().begin(); y != r.cols().end(); ++y) {
128 double dist = distance_squared(s, my_map[x][y]);
129 if (dist < min_distance_squared) {
130 min_distance_squared = dist;
131 min_x = x;
132 min_y = y;
133 }
134 if (cancel_test && oneapi::tbb::is_current_task_group_canceling()) {
135 xval = r.rows().begin();
136 yval = r.cols().begin();
137 return DBL_MAX;
138 }
139 }
140 }
141 xval = min_x;
142 yval = min_y;
143 return sqrt(min_distance_squared);
144 }
145
epoch_update_range(SOM_element const & s,int epoch,int min_x,int min_y,double radius,double learning_rate,oneapi::tbb::blocked_range<int> & r)146 void SOMap::epoch_update_range(SOM_element const &s,
147 int epoch,
148 int min_x,
149 int min_y,
150 double radius,
151 double learning_rate,
152 oneapi::tbb::blocked_range<int> &r) {
153 int min_xiter = (int)((double)min_x - radius);
154 if (min_xiter < 0)
155 min_xiter = 0;
156 int max_xiter = (int)((double)min_x + radius);
157 if (max_xiter > (int)my_map.size() - 1)
158 max_xiter = (int)my_map.size() - 1;
159 for (int xx = r.begin(); xx <= r.end(); ++xx) {
160 double xrsq = (xx - min_x) * (xx - min_x);
161 double ysq = radius * radius - xrsq; // max extent of y influence
162 double yd;
163 if (ysq > 0) {
164 yd = sqrt(ysq);
165 int lb = (int)(min_y - yd);
166 int ub = (int)(min_y + yd);
167 for (int yy = lb; yy < ub; ++yy) {
168 if (yy >= 0 && yy < (int)my_map[xx].size()) {
169 // [xx, yy] is in the range of the update.
170 double my_rsq = xrsq + (yy - min_y) * (yy - min_y); // distance from BMU squared
171 double theta = exp(-(radius * radius) / (2.0 * my_rsq));
172 add_fraction_of_difference(my_map[xx][yy], s, theta * learning_rate);
173 }
174 }
175 }
176 }
177 }
178
teach(teaching_vector_type & in)179 void SOMap::teach(teaching_vector_type &in) {
180 for (int i = 0; i < nPasses; ++i) {
181 int j = (int)(randval(0, (double)in.size())); // this won't be reproducible.
182 if (j == in.size())
183 --j;
184
185 int min_x = -1;
186 int min_y = -1;
187 subsquare_type br2(0, (int)my_map.size(), 1, 0, (int)my_map[0].size(), 1);
188 (void)BMU_range(in[j], min_x, min_y, br2); // just need min_x, min_y
189 // radius of interest
190 double radius = max_radius * exp(-(double)i * radius_decay_rate);
191 // update circle is min_xiter to max_xiter inclusive.
192 double learning_rate = max_learning_rate * exp(-(double)i * learning_decay_rate);
193 epoch_update(in[j], i, min_x, min_y, radius, learning_rate);
194 }
195 }
196
debug_output()197 void SOMap::debug_output() {
198 printf("SOMap:\n");
199 for (int i = 0; i < (int)(this->my_map.size()); ++i) {
200 for (int j = 0; j < (int)(this->my_map[i].size()); ++j) {
201 printf("map[%d, %d] == ", i, j);
202 remark_SOM_element(this->my_map[i][j]);
203 printf("\n");
204 }
205 }
206 }
207
208 #define RED 0
209 #define GREEN 1
210 #define BLUE 2
211
readInputData()212 void readInputData() {
213 my_teaching.push_back(SOM_element());
214 my_teaching.push_back(SOM_element());
215 my_teaching.push_back(SOM_element());
216 my_teaching.push_back(SOM_element());
217 my_teaching.push_back(SOM_element());
218 my_teaching[0][RED] = 1.0;
219 my_teaching[0][GREEN] = 0.0;
220 my_teaching[0][BLUE] = 0.0;
221 my_teaching[1][RED] = 0.0;
222 my_teaching[1][GREEN] = 1.0;
223 my_teaching[1][BLUE] = 0.0;
224 my_teaching[2][RED] = 0.0;
225 my_teaching[2][GREEN] = 0.0;
226 my_teaching[2][BLUE] = 1.0;
227 my_teaching[3][RED] = 0.3;
228 my_teaching[3][GREEN] = 0.3;
229 my_teaching[3][BLUE] = 0.0;
230 my_teaching[4][RED] = 0.5;
231 my_teaching[4][GREEN] = 0.5;
232 my_teaching[4][BLUE] = 0.9;
233 }
234