kohonen_som_topology.cpp

Go to the documentation of this file.

 
#define _USE_MATH_DEFINES  //< required for MS Visual C++
#include <algorithm>
#include <array>
#include <cerrno>
#include <cmath>
#include <cstdlib>
#include <cstring>
#include <ctime>
#include <fstream>
#include <iostream>
#include <valarray>
#include <vector>
#ifdef _OPENMP  // check if OpenMP based parallellization is available
#include <omp.h>
#endif
 
double _random(double a, double b) {
    return ((b - a) * (std::rand() % 100) / 100.f) + a;
}
 
int save_2d_data(const char *fname,
                 const std::vector<std::valarray<double>> &X) {
    size_t num_points = X.size();       // number of rows
    size_t num_features = X[0].size();  // number of columns
 
    std::ofstream fp;
    fp.open(fname);
    if (!fp.is_open()) {
        // error with opening file to write
        std::cerr << "Error opening file " << fname << ": "
                  << std::strerror(errno) << "\n";
        return -1;
    }
 
    // for each point in the array
    for (int i = 0; i < num_points; i++) {
        // for each feature in the array
        for (int j = 0; j < num_features; j++) {
            fp << X[i][j];               // print the feature value
            if (j < num_features - 1) {  // if not the last feature
                fp << ",";               // suffix comma
            }
        }
        if (i < num_points - 1) {  // if not the last row
            fp << "\n";            // start a new line
        }
    }
 
    fp.close();
    return 0;
}
 
void get_min_2d(const std::vector<std::valarray<double>> &X, double *val,
                int *x_idx, int *y_idx) {
    val[0] = INFINITY;  // initial min value
    size_t N = X.size();
 
    for (int i = 0; i < N; i++) {  // traverse each x-index
        auto result = std::min_element(std::begin(X[i]), std::end(X[i]));
        double d_min = *result;
        std::ptrdiff_t j = std::distance(std::begin(X[i]), result);
 
        if (d_min < val[0]) {  // if a lower value is found
                               // save the value and its index
            x_idx[0] = i;
            y_idx[0] = j;
            val[0] = d_min;
        }
    }
}
 
namespace machine_learning {
constexpr double MIN_DISTANCE = 1e-4;
 
int save_u_matrix(const char *fname,
                  const std::vector<std::vector<std::valarray<double>>> &W) {
    std::ofstream fp(fname);
    if (!fp) {  // error with fopen
        std::cerr << "File error (" << fname << "): " << std::strerror(errno)
                  << std::endl;
        return -1;
    }
 
    // neighborhood range
    unsigned int R = 1;
 
    for (int i = 0; i < W.size(); i++) {         // for each x
        for (int j = 0; j < W[0].size(); j++) {  // for each y
            double distance = 0.f;
 
            int from_x = std::max<int>(0, i - R);
            int to_x = std::min<int>(W.size(), i + R + 1);
            int from_y = std::max<int>(0, j - R);
            int to_y = std::min<int>(W[0].size(), j + R + 1);
            int l = 0, m = 0;
#ifdef _OPENMP
#pragma omp parallel for reduction(+ : distance)
#endif
            for (l = from_x; l < to_x; l++) {      // scan neighborhoor in x
                for (m = from_y; m < to_y; m++) {  // scan neighborhood in y
                    auto d = W[i][j] - W[l][m];
                    double d2 = std::pow(d, 2).sum();
                    distance += std::sqrt(d2);
                    // distance += d2;
                }
            }
 
            distance /= R * R;          // mean distance from neighbors
            fp << distance;             // print the mean separation
            if (j < W[0].size() - 1) {  // if not the last column
                fp << ',';              // suffix comma
            }
        }
        if (i < W.size() - 1) {  // if not the last row
            fp << '\n';          // start a new line
        }
    }
 
    fp.close();
    return 0;
}
 
double update_weights(const std::valarray<double> &X,
                      std::vector<std::vector<std::valarray<double>>> *W,
                      std::vector<std::valarray<double>> *D, double alpha,
                      int R) {
    int x = 0, y = 0;
    int num_out_x = static_cast<int>(W->size());       // output nodes - in X
    int num_out_y = static_cast<int>(W[0][0].size());  // output nodes - in Y
    // int num_features = static_cast<int>(W[0][0][0].size());  //  features =
    // in Z
    double d_min = 0.f;
 
#ifdef _OPENMP
#pragma omp for
#endif
    // step 1: for each output point
    for (x = 0; x < num_out_x; x++) {
        for (y = 0; y < num_out_y; y++) {
            (*D)[x][y] = 0.f;
            // compute Euclidian distance of each output
            // point from the current sample
            auto d = ((*W)[x][y] - X);
            (*D)[x][y] = (d * d).sum();
            (*D)[x][y] = std::sqrt((*D)[x][y]);
        }
    }
 
    // step 2:  get closest node i.e., node with snallest Euclidian distance
    // to the current pattern
    int d_min_x = 0, d_min_y = 0;
    get_min_2d(*D, &d_min, &d_min_x, &d_min_y);
 
    // step 3a: get the neighborhood range
    int from_x = std::max(0, d_min_x - R);
    int to_x = std::min(num_out_x, d_min_x + R + 1);
    int from_y = std::max(0, d_min_y - R);
    int to_y = std::min(num_out_y, d_min_y + R + 1);
 
    // step 3b: update the weights of nodes in the
    // neighborhood
#ifdef _OPENMP
#pragma omp for
#endif
    for (x = from_x; x < to_x; x++) {
        for (y = from_y; y < to_y; y++) {
            /* you can enable the following normalization if needed.
   personally, I found it detrimental to convergence */
            // const double s2pi = sqrt(2.f * M_PI);
            // double normalize = 1.f / (alpha * s2pi);
 
            /* apply scaling inversely proportional to distance from the
               current node */
            double d2 =
                (d_min_x - x) * (d_min_x - x) + (d_min_y - y) * (d_min_y - y);
            double scale_factor = std::exp(-d2 / (2.f * alpha * alpha));
 
            (*W)[x][y] += (X - (*W)[x][y]) * alpha * scale_factor;
        }
    }
    return d_min;
}
 
void kohonen_som(const std::vector<std::valarray<double>> &X,
                 std::vector<std::vector<std::valarray<double>>> *W,
                 double alpha_min) {
    size_t num_samples = X.size();  // number of rows
    // size_t num_features = X[0].size();  // number of columns
    size_t num_out = W->size();  // output matrix size
    size_t R = num_out >> 2, iter = 0;
    double alpha = 1.f;
 
    std::vector<std::valarray<double>> D(num_out);
    for (int i = 0; i < num_out; i++) D[i] = std::valarray<double>(num_out);
 
    double dmin = 1.f;        // average minimum distance of all samples
    double past_dmin = 1.f;   // average minimum distance of all samples
    double dmin_ratio = 1.f;  // change per step
 
    // Loop alpha from 1 to slpha_min
    for (; alpha > 0 && dmin_ratio > 1e-5; alpha -= 1e-4, iter++) {
        // Loop for each sample pattern in the data set
        for (int sample = 0; sample < num_samples; sample++) {
            // update weights for the current input pattern sample
            dmin += update_weights(X[sample], W, &D, alpha, R);
        }
 
        // every 100th iteration, reduce the neighborhood range
        if (iter % 300 == 0 && R > 1) {
            R--;
        }
 
        dmin /= num_samples;
 
        // termination condition variable -> % change in minimum distance
        dmin_ratio = (past_dmin - dmin) / past_dmin;
        if (dmin_ratio < 0) {
            dmin_ratio = 1.f;
        }
        past_dmin = dmin;
 
        std::cout << "iter: " << iter << "\t alpha: " << alpha << "\t R: " << R
                  << "\t d_min: " << dmin_ratio << "\r";
    }
 
    std::cout << "\n";
}
 
}  // namespace machine_learning
 
using machine_learning::kohonen_som;
using machine_learning::save_u_matrix;
 
void test_2d_classes(std::vector<std::valarray<double>> *data) {
    const int N = data->size();
    const double R = 0.3;  // radius of cluster
    int i = 0;
    const int num_classes = 4;
    std::array<std::array<double, 2>, num_classes> centres = {
        // centres of each class cluster
        std::array<double, 2>({.5, .5}),   // centre of class 1
        std::array<double, 2>({.5, -.5}),  // centre of class 2
        std::array<double, 2>({-.5, .5}),  // centre of class 3
        std::array<double, 2>({-.5, -.5})  // centre of class 4
    };
 
#ifdef _OPENMP
#pragma omp for
#endif
    for (i = 0; i < N; i++) {
        // select a random class for the point
        int cls = std::rand() % num_classes;
 
        // create random coordinates (x,y,z) around the centre of the class
        data[0][i][0] = _random(centres[cls][0] - R, centres[cls][0] + R);
        data[0][i][1] = _random(centres[cls][1] - R, centres[cls][1] + R);
 
        /* The follosing can also be used
        for (int j = 0; j < 2; j++)
            data[i][j] = _random(centres[class][j] - R, centres[class][j] + R);
        */
    }
}
 
void test1() {
    int j = 0, N = 300;
    int features = 2;
    int num_out = 30;
    std::vector<std::valarray<double>> X(N);
    std::vector<std::vector<std::valarray<double>>> W(num_out);
    for (int i = 0; i < std::max(num_out, N); i++) {
        // loop till max(N, num_out)
        if (i < N) {  // only add new arrays if i < N
            X[i] = std::valarray<double>(features);
        }
        if (i < num_out) {  // only add new arrays if i < num_out
            W[i] = std::vector<std::valarray<double>>(num_out);
            for (int k = 0; k < num_out; k++) {
                W[i][k] = std::valarray<double>(features);
#ifdef _OPENMP
#pragma omp for
#endif
                for (j = 0; j < features; j++) {
                    // preallocate with random initial weights
                    W[i][k][j] = _random(-10, 10);
                }
            }
        }
    }
 
    test_2d_classes(&X);  // create test data around circumference of a circle
    save_2d_data("test1.csv", X);  // save test data points
    save_u_matrix("w11.csv", W);   // save initial random weights
    kohonen_som(X, &W, 1e-4);      // train the SOM
    save_u_matrix("w12.csv", W);   // save the resultant weights
}
 
void test_3d_classes1(std::vector<std::valarray<double>> *data) {
    const size_t N = data->size();
    const double R = 0.3;  // radius of cluster
    int i = 0;
    const int num_classes = 4;
    const std::array<std::array<double, 3>, num_classes> centres = {
        // centres of each class cluster
        std::array<double, 3>({.5, .5, .5}),    // centre of class 1
        std::array<double, 3>({.5, -.5, -.5}),  // centre of class 2
        std::array<double, 3>({-.5, .5, .5}),   // centre of class 3
        std::array<double, 3>({-.5, -.5 - .5})  // centre of class 4
    };
 
#ifdef _OPENMP
#pragma omp for
#endif
    for (i = 0; i < N; i++) {
        // select a random class for the point
        int cls = std::rand() % num_classes;
 
        // create random coordinates (x,y,z) around the centre of the class
        data[0][i][0] = _random(centres[cls][0] - R, centres[cls][0] + R);
        data[0][i][1] = _random(centres[cls][1] - R, centres[cls][1] + R);
        data[0][i][2] = _random(centres[cls][2] - R, centres[cls][2] + R);
 
        /* The follosing can also be used
        for (int j = 0; j < 3; j++)
            data[i][j] = _random(centres[class][j] - R, centres[class][j] + R);
        */
    }
}
 
void test2() {
    int j = 0, N = 300;
    int features = 3;
    int num_out = 30;
    std::vector<std::valarray<double>> X(N);
    std::vector<std::vector<std::valarray<double>>> W(num_out);
    for (int i = 0; i < std::max(num_out, N); i++) {
        // loop till max(N, num_out)
        if (i < N) {  // only add new arrays if i < N
            X[i] = std::valarray<double>(features);
        }
        if (i < num_out) {  // only add new arrays if i < num_out
            W[i] = std::vector<std::valarray<double>>(num_out);
            for (int k = 0; k < num_out; k++) {
                W[i][k] = std::valarray<double>(features);
#ifdef _OPENMP
#pragma omp for
#endif
                for (j = 0; j < features; j++) {
                    // preallocate with random initial weights
                    W[i][k][j] = _random(-10, 10);
                }
            }
        }
    }
 
    test_3d_classes1(&X);  // create test data around circumference of a circle
    save_2d_data("test2.csv", X);  // save test data points
    save_u_matrix("w21.csv", W);   // save initial random weights
    kohonen_som(X, &W, 1e-4);      // train the SOM
    save_u_matrix("w22.csv", W);   // save the resultant weights
}
 
void test_3d_classes2(std::vector<std::valarray<double>> *data) {
    const size_t N = data->size();
    const double R = 0.2;  // radius of cluster
    int i = 0;
    const int num_classes = 8;
    const std::array<std::array<double, 3>, num_classes> centres = {
        // centres of each class cluster
        std::array<double, 3>({.5, .5, .5}),    // centre of class 1
        std::array<double, 3>({.5, .5, -.5}),   // centre of class 2
        std::array<double, 3>({.5, -.5, .5}),   // centre of class 3
        std::array<double, 3>({.5, -.5, -.5}),  // centre of class 4
        std::array<double, 3>({-.5, .5, .5}),   // centre of class 5
        std::array<double, 3>({-.5, .5, -.5}),  // centre of class 6
        std::array<double, 3>({-.5, -.5, .5}),  // centre of class 7
        std::array<double, 3>({-.5, -.5, -.5})  // centre of class 8
    };
 
#ifdef _OPENMP
#pragma omp for
#endif
    for (i = 0; i < N; i++) {
        // select a random class for the point
        int cls = std::rand() % num_classes;
 
        // create random coordinates (x,y,z) around the centre of the class
        data[0][i][0] = _random(centres[cls][0] - R, centres[cls][0] + R);
        data[0][i][1] = _random(centres[cls][1] - R, centres[cls][1] + R);
        data[0][i][2] = _random(centres[cls][2] - R, centres[cls][2] + R);
 
        /* The follosing can also be used
        for (int j = 0; j < 3; j++)
            data[i][j] = _random(centres[class][j] - R, centres[class][j] + R);
        */
    }
}
 
void test3() {
    int j = 0, N = 500;
    int features = 3;
    int num_out = 30;
    std::vector<std::valarray<double>> X(N);
    std::vector<std::vector<std::valarray<double>>> W(num_out);
    for (int i = 0; i < std::max(num_out, N); i++) {
        // loop till max(N, num_out)
        if (i < N) {  // only add new arrays if i < N
            X[i] = std::valarray<double>(features);
        }
        if (i < num_out) {  // only add new arrays if i < num_out
            W[i] = std::vector<std::valarray<double>>(num_out);
            for (int k = 0; k < num_out; k++) {
                W[i][k] = std::valarray<double>(features);
#ifdef _OPENMP
#pragma omp for
#endif
                for (j = 0; j < features; j++) {
                    // preallocate with random initial weights
                    W[i][k][j] = _random(-10, 10);
                }
            }
        }
    }
 
    test_3d_classes2(&X);  // create test data around circumference of a circle
    save_2d_data("test3.csv", X);  // save test data points
    save_u_matrix("w31.csv", W);   // save initial random weights
    kohonen_som(X, &W, 1e-4);      // train the SOM
    save_u_matrix("w32.csv", W);   // save the resultant weights
}
 
double get_clock_diff(clock_t start_t, clock_t end_t) {
    return static_cast<double>(end_t - start_t) / CLOCKS_PER_SEC;
}
 
int main(int argc, char **argv) {
#ifdef _OPENMP
    std::cout << "Using OpenMP based parallelization\n";
#else
    std::cout << "NOT using OpenMP based parallelization\n";
#endif
 
    std::srand(std::time(nullptr));
 
    std::clock_t start_clk = std::clock();
    test1();
    auto end_clk = std::clock();
    std::cout << "Test 1 completed in " << get_clock_diff(start_clk, end_clk)
              << " sec\n";
 
    start_clk = std::clock();
    test2();
    end_clk = std::clock();
    std::cout << "Test 2 completed in " << get_clock_diff(start_clk, end_clk)
              << " sec\n";
 
    start_clk = std::clock();
    test3();
    end_clk = std::clock();
    std::cout << "Test 3 completed in " << get_clock_diff(start_clk, end_clk)
              << " sec\n";
 
    std::cout
        << "(Note: Calculated times include: creating test sets, training "
           "model and writing files to disk.)\n\n";
    return 0;
}