kohonen_som_trace.cpp

Go to the documentation of this file.

 
#define _USE_MATH_DEFINES  // required for MS Visual C++
#include <algorithm>
#include <array>
#include <cmath>
#include <cstdlib>
#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_nd_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 << "\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;
}
 
namespace machine_learning {
 
void update_weights(const std::valarray<double> &x,
                    std::vector<std::valarray<double>> *W,
                    std::valarray<double> *D, double alpha, int R) {
    int j = 0, k = 0;
    int num_out = W->size();  // number of SOM output nodes
    // int num_features = x.size();  // number of data features
 
#ifdef _OPENMP
#pragma omp for
#endif
    // step 1: for each output point
    for (j = 0; j < num_out; j++) {
        // compute Euclidian distance of each output
        // point from the current sample
        (*D)[j] = (((*W)[j] - x) * ((*W)[j] - x)).sum();
    }
 
    // step 2:  get closest node i.e., node with snallest Euclidian distance to
    // the current pattern
    auto result = std::min_element(std::begin(*D), std::end(*D));
    // double d_min = *result;
    int d_min_idx = std::distance(std::begin(*D), result);
 
    // step 3a: get the neighborhood range
    int from_node = std::max(0, d_min_idx - R);
    int to_node = std::min(num_out, d_min_idx + R + 1);
 
    // step 3b: update the weights of nodes in the
    // neighborhood
#ifdef _OPENMP
#pragma omp for
#endif
    for (j = from_node; j < to_node; j++) {
        // update weights of nodes in the neighborhood
        (*W)[j] += alpha * (x - (*W)[j]);
    }
}
 
void kohonen_som_tracer(const std::vector<std::valarray<double>> &X,
                        std::vector<std::valarray<double>> *W,
                        double alpha_min) {
    int num_samples = X.size();  // number of rows
    // int num_features = X[0].size();  // number of columns
    int num_out = W->size();  // number of rows
    int R = num_out >> 2, iter = 0;
    double alpha = 1.f;
 
    std::valarray<double> D(num_out);
 
    // Loop alpha from 1 to slpha_min
    do {
        // 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
            update_weights(X[sample], W, &D, alpha, R);
        }
 
        // every 10th iteration, reduce the neighborhood range
        if (iter % 10 == 0 && R > 1) {
            R--;
        }
 
        alpha -= 0.01;
        iter++;
    } while (alpha > alpha_min);
}
 
}  // namespace machine_learning
 
using machine_learning::kohonen_som_tracer;
 
void test_circle(std::vector<std::valarray<double>> *data) {
    const int N = data->size();
    const double R = 0.75, dr = 0.3;
    double a_t = 0., b_t = 2.f * M_PI;  // theta random between 0 and 2*pi
    double a_r = R - dr, b_r = R + dr;  // radius random between R-dr and R+dr
    int i = 0;
 
#ifdef _OPENMP
#pragma omp for
#endif
    for (i = 0; i < N; i++) {
        double r = _random(a_r, b_r);      // random radius
        double theta = _random(a_t, b_t);  // random theta
        data[0][i][0] = r * cos(theta);    // convert from polar to cartesian
        data[0][i][1] = r * sin(theta);
    }
}
 
void test1() {
    int j = 0, N = 500;
    int features = 2;
    int num_out = 50;
    std::vector<std::valarray<double>> X(N);
    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::valarray<double>(features);
 
#ifdef _OPENMP
#pragma omp for
#endif
            for (j = 0; j < features; j++) {
                // preallocate with random initial weights
                W[i][j] = _random(-1, 1);
            }
        }
    }
 
    test_circle(&X);  // create test data around circumference of a circle
    save_nd_data("test1.csv", X);    // save test data points
    save_nd_data("w11.csv", W);      // save initial random weights
    kohonen_som_tracer(X, &W, 0.1);  // train the SOM
    save_nd_data("w12.csv", W);      // save the resultant weights
}
 
void test_lamniscate(std::vector<std::valarray<double>> *data) {
    const int N = data->size();
    const double dr = 0.2;
    int i = 0;
 
#ifdef _OPENMP
#pragma omp for
#endif
    for (i = 0; i < N; i++) {
        double dx = _random(-dr, dr);     // random change in x
        double dy = _random(-dr, dr);     // random change in y
        double theta = _random(0, M_PI);  // random theta
        data[0][i][0] = dx + cos(theta);  // convert from polar to cartesian
        data[0][i][1] = dy + sin(2. * theta) / 2.f;
    }
}
 
void test2() {
    int j = 0, N = 500;
    int features = 2;
    int num_out = 20;
    std::vector<std::valarray<double>> X(N);
    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::valarray<double>(features);
 
#ifdef _OPENMP
#pragma omp for
#endif
            for (j = 0; j < features; j++) {
                // preallocate with random initial weights
                W[i][j] = _random(-1, 1);
            }
        }
    }
 
    test_lamniscate(&X);              // create test data around the lamniscate
    save_nd_data("test2.csv", X);     // save test data points
    save_nd_data("w21.csv", W);       // save initial random weights
    kohonen_som_tracer(X, &W, 0.01);  // train the SOM
    save_nd_data("w22.csv", W);       // save the resultant weights
}
 
void test_3d_classes(std::vector<std::valarray<double>> *data) {
    const int N = data->size();
    const double R = 0.1;  // radius of cluster
    int i = 0;
    const int num_classes = 8;
    const std::array<const std::array<double, 3>, num_classes> centres = {
        // centres of each class cluster
        std::array<double, 3>({.5, .5, .5}),    // centre of class 0
        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
    };
 
#ifdef _OPENMP
#pragma omp for
#endif
    for (i = 0; i < N; i++) {
        int cls =
            std::rand() % num_classes;  // select a random class for the point
 
        // 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[0][i][j] = _random(centres[cls][j] - R, centres[cls][j] + R);
        */
    }
}
 
void test3() {
    int j = 0, N = 200;
    int features = 3;
    int num_out = 20;
    std::vector<std::valarray<double>> X(N);
    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::valarray<double>(features);
 
#ifdef _OPENMP
#pragma omp for
#endif
            for (j = 0; j < features; j++) {
                // preallocate with random initial weights
                W[i][j] = _random(-1, 1);
            }
        }
    }
 
    test_3d_classes(&X);              // create test data around the lamniscate
    save_nd_data("test3.csv", X);     // save test data points
    save_nd_data("w31.csv", W);       // save initial random weights
    kohonen_som_tracer(X, &W, 0.01);  // train the SOM
    save_nd_data("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;
}