examples/graph/som/som.cpp

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