diff --git a/DIRECTORY.md b/DIRECTORY.md
index 35de16ae..39549c55 100644
--- a/DIRECTORY.md
+++ b/DIRECTORY.md
@@ -203,6 +203,7 @@
 
 ## Machine Learning
   * [Adaline Learning](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/machine_learning/adaline_learning.c)
+  * [Kohonen Som](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/machine_learning/kohonen_som.c)
 
 ## Misc
   * [Armstrong Number](https://github.com/TheAlgorithms/C-Plus-Plus/blob/master/misc/armstrong_number.c)
diff --git a/machine_learning/adaline_learning.c b/machine_learning/adaline_learning.c
index 7c69abef..3bac6823 100644
--- a/machine_learning/adaline_learning.c
+++ b/machine_learning/adaline_learning.c
@@ -168,7 +168,7 @@ double fit_sample(struct adaline *ada, const double *x, const int y)
  * \param[in] y known output value for each feature vector
  * \param[in] N number of training samples
  */
-void fit(struct adaline *ada, const double **X, const int *y, const int N)
+void fit(struct adaline *ada, double **X, const int *y, const int N)
 {
     double avg_pred_error = 1.f;
 
diff --git a/machine_learning/kohonen_som.c b/machine_learning/kohonen_som.c
new file mode 100644
index 00000000..d0415e9f
--- /dev/null
+++ b/machine_learning/kohonen_som.c
@@ -0,0 +1,486 @@
+/**
+ * \file
+ * \brief [Kohonen self organizing
+ * map](https://en.wikipedia.org/wiki/Self-organizing_map) (1D)
+ *
+ * This example implements a powerful self organizing map algorithm in 1D.
+ * The algorithm creates a connected network of weights that closely
+ * follows the given data points. This this creates a chain of nodes that
+ * resembles the given input shape.
+ */
+#define _USE_MATH_DEFINES // required for MS Visual C
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <time.h>
+#ifdef _OPENMP // check if OpenMP based parallellization is available
+#include <omp.h>
+#endif
+
+/**
+ * Helper function to generate a random number in a given interval.
+ * \n Steps:
+ * 1. `r1 = rand() % 100` gets a random number between 0 and 99
+ * 2. `r2 = r1 / 100` converts random number to be between 0 and 0.99
+ * 3. scale and offset the random number to given range of \f$[a,b]\f$
+ *
+ * \param[in] a lower limit
+ * \param[in] b upper limit
+ * \returns random number in the range \f$[a,b]\f$
+ */
+double _random(double a, double b)
+{
+    return ((b - a) * (rand() % 100) / 100.f) + a;
+}
+
+/**
+ * Save a given n-dimensional data martix to file.
+ *
+ * \param[in] fname filename to save in (gets overwriten without confirmation)
+ * \param[in] X matrix to save
+ * \param[in] num_points rows in the matrix = number of points
+ * \param[in] num_features columns in the matrix = dimensions of points
+ * \returns 0 if all ok
+ * \returns -1 if file creation failed
+ */
+int save_nd_data(const char *fname, double **X, int num_points,
+                 int num_features)
+{
+    FILE *fp = fopen(fname, "wt");
+    if (!fp) // error with fopen
+    {
+        char msg[120];
+        sprintf(msg, "File error (%s): ", fname);
+        perror(msg);
+        return -1;
+    }
+
+    for (int i = 0; i < num_points; i++) // for each point in the array
+    {
+        for (int j = 0; j < num_features; j++) // for each feature in the array
+        {
+            fprintf(fp, "%.4g", X[i][j]); // print the feature value
+            if (j < num_features - 1)     // if not the last feature
+                fprintf(fp, ",");         // suffix comma
+        }
+        if (i < num_points - 1) // if not the last row
+            fprintf(fp, "\n");  // start a new line
+    }
+    fclose(fp);
+    return 0;
+}
+
+/**
+ * Get minimum value and index of the value in a vector
+ * \param[in] x vector to search
+ * \param[in] N number of points in the vector
+ * \param[out] val minimum value found
+ * \param[out] idx index where minimum value was found
+ */
+void get_min_1d(double const *X, int N, double *val, int *idx)
+{
+    val[0] = INFINITY; // initial min value
+
+    for (int i = 0; i < N; i++) // check each value
+    {
+        if (X[i] < val[0]) // if a lower value is found
+        {                  // save the value and its index
+            idx[0] = i;
+            val[0] = X[i];
+        }
+    }
+}
+
+/**
+ * Update weights of the SOM using Kohonen algorithm
+ *
+ * \param[in] X data point
+ * \param[in,out] W weights matrix
+ * \param[in,out] D temporary vector to store distances
+ * \param[in] num_out number of output points
+ * \param[in] num_features number of features per input sample
+ * \param[in] alpha learning rate \f$0<\alpha\le1\f$
+ * \param[in] R neighborhood range
+ */
+void update_weights(double const *x, double *const *W, double *D, int num_out,
+                    int num_features, double alpha, int R)
+{
+    int j, k;
+
+#ifdef _OPENMP
+#pragma omp for
+#endif
+    // step 1: for each output point
+    for (j = 0; j < num_out; j++)
+    {
+        D[j] = 0.f;
+        // compute Euclidian distance of each output
+        // point from the current sample
+        for (k = 0; k < num_features; k++)
+            D[j] += (W[j][k] - x[k]) * (W[j][k] - x[k]);
+    }
+
+    // step 2:  get closest node i.e., node with snallest Euclidian distance to
+    // the current pattern
+    int d_min_idx;
+    double d_min;
+    get_min_1d(D, num_out, &d_min, &d_min_idx);
+
+    // step 3a: get the neighborhood range
+    int from_node = 0 > (d_min_idx - R) ? 0 : d_min_idx - R;
+    int to_node = num_out < (d_min_idx + R + 1) ? 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++)
+        for (k = 0; k < num_features; k++)
+            // update weights of nodes in the neighborhood
+            W[j][k] += alpha * (x[k] - W[j][k]);
+}
+
+/**
+ * Apply incremental algorithm with updating neighborhood and learning rates
+ * on all samples in the given datset.
+ *
+ * \param[in] X data set
+ * \param[in,out] W weights matrix
+ * \param[in] D temporary vector to store distances
+ * \param[in] num_samples number of output points
+ * \param[in] num_features number of features per input sample
+ * \param[in] num_out number of output points
+ * \param[in] alpha_min terminal value of alpha
+ */
+void kohonen_som_tracer(double **X, double *const *W, int num_samples,
+                        int num_features, int num_out, double alpha_min)
+{
+    int R = num_out >> 2, iter = 0;
+    double alpha = 1.f;
+    double *D = (double *)malloc(num_out * sizeof(double));
+
+    // Loop alpha from 1 to slpha_min
+    for (; alpha > alpha_min; alpha -= 0.01, iter++)
+    {
+        // Loop for each sample pattern in the data set
+        for (int sample = 0; sample < num_samples; sample++)
+        {
+            const double *x = X[sample];
+            // update weights for the current input pattern sample
+            update_weights(x, W, D, num_out, num_features, alpha, R);
+        }
+
+        // every 10th iteration, reduce the neighborhood range
+        if (iter % 10 == 0 && R > 1)
+            R--;
+    }
+
+    free(D);
+}
+
+/** Creates a random set of points distributed *near* the circumference
+ * of a circle and trains an SOM that finds that circular pattern. The
+ * generating function is
+ * \f{eqnarray*}{ \f}
+ *
+ * \param[out] data matrix to store data in
+ * \param[in] N number of points required
+ */
+void test_circle(double *const *data, int N)
+{
+    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;
+
+#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[i][0] = r * cos(theta);      // convert from polar to cartesian
+        data[i][1] = r * sin(theta);
+    }
+}
+
+/** Test that creates a random set of points distributed *near* the
+ * circumference of a circle and trains an SOM that finds that circular pattern.
+ * The following [CSV](https://en.wikipedia.org/wiki/Comma-separated_values)
+ * files are created to validate the execution:
+ * * `test1.csv`: random test samples points with a circular pattern
+ * * `w11.csv`: initial random map
+ * * `w12.csv`: trained SOM map
+ *
+ * The outputs can be readily plotted in [gnuplot](https:://gnuplot.info) using
+ * the following snippet
+ * ```gnuplot
+ * set datafile separator ','
+ * plot "test1.csv" title "original", \
+ *      "w11.csv" title "w1", \
+ *      "w12.csv" title "w2"
+ * ```
+ */
+void test1()
+{
+    int j, N = 500;
+    int features = 2;
+    int num_out = 50;
+    double **X = (double **)malloc(N * sizeof(double *));
+    double **W = (double **)malloc(num_out * sizeof(double *));
+    for (int i = 0; i < (num_out > N ? num_out : N);
+         i++) // loop till max(N, num_out)
+    {
+        if (i < N) // only add new arrays if i < N
+            X[i] = (double *)malloc(features * sizeof(double));
+        if (i < num_out) // only add new arrays if i < num_out
+        {
+            W[i] = (double *)malloc(features * sizeof(double));
+#ifdef _OPENMP
+#pragma omp for
+#endif
+            // preallocate with random initial weights
+            for (j = 0; j < features; j++)
+                W[i][j] = _random(-1, 1);
+        }
+    }
+
+    test_circle(X, N); // create test data around circumference of a circle
+    save_nd_data("test1.csv", X, N, features); // save test data points
+    save_nd_data("w11.csv", W, num_out,
+                 features); // save initial random weights
+    kohonen_som_tracer(X, W, N, features, num_out, 0.1); // train the SOM
+    save_nd_data("w12.csv", W, num_out, features); // save the resultant weights
+
+    for (int i = 0; i < (num_out > N ? num_out : N); i++)
+    {
+        if (i < N)
+            free(X[i]);
+        if (i < num_out)
+            free(W[i]);
+    }
+}
+
+/** Creates a random set of points distributed *near* the locus
+ * of the [Lamniscate of
+ * Gerono](https://en.wikipedia.org/wiki/Lemniscate_of_Gerono) and trains an SOM
+ * that finds that circular pattern. \param[out] data matrix to store data in
+ * \param[in] N number of points required
+ */
+void test_lamniscate(double *const *data, int N)
+{
+    const double dr = 0.2;
+    int i;
+
+#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[i][0] = dx + cos(theta);    // convert from polar to cartesian
+        data[i][1] = dy + sin(2. * theta) / 2.f;
+    }
+}
+
+/** Test that creates a random set of points distributed *near* the locus
+ * of the [Lamniscate of
+ * Gerono](https://en.wikipedia.org/wiki/Lemniscate_of_Gerono) and trains an SOM
+ * that finds that circular pattern. The following
+ * [CSV](https://en.wikipedia.org/wiki/Comma-separated_values) files are created
+ * to validate the execution:
+ * * `test2.csv`: random test samples points with a lamniscate pattern
+ * * `w21.csv`: initial random map
+ * * `w22.csv`: trained SOM map
+ *
+ * The outputs can be readily plotted in [gnuplot](https:://gnuplot.info) using
+ * the following snippet
+ * ```gnuplot
+ * set datafile separator ','
+ * plot "test2.csv" title "original", \
+ *      "w21.csv" title "w1", \
+ *      "w22.csv" title "w2"
+ * ```
+ */
+void test2()
+{
+    int j, N = 500;
+    int features = 2;
+    int num_out = 20;
+    double **X = (double **)malloc(N * sizeof(double *));
+    double **W = (double **)malloc(num_out * sizeof(double *));
+    for (int i = 0; i < (num_out > N ? num_out : N); i++)
+    {
+        if (i < N) // only add new arrays if i < N
+            X[i] = (double *)malloc(features * sizeof(double));
+        if (i < num_out) // only add new arrays if i < num_out
+        {
+            W[i] = (double *)malloc(features * sizeof(double));
+
+#ifdef _OPENMP
+#pragma omp for
+#endif
+            // preallocate with random initial weights
+            for (j = 0; j < features; j++)
+                W[i][j] = _random(-1, 1);
+        }
+    }
+
+    test_lamniscate(X, N); // create test data around the lamniscate
+    save_nd_data("test2.csv", X, N, features); // save test data points
+    save_nd_data("w21.csv", W, num_out,
+                 features); // save initial random weights
+    kohonen_som_tracer(X, W, N, features, num_out, 0.01); // train the SOM
+    save_nd_data("w22.csv", W, num_out, features); // save the resultant weights
+
+    for (int i = 0; i < (num_out > N ? num_out : N); i++)
+    {
+        if (i < N)
+            free(X[i]);
+        if (i < num_out)
+            free(W[i]);
+    }
+    free(X);
+    free(W);
+}
+
+/** Creates a random set of points distributed *near* the locus
+ * of the [Lamniscate of
+ * Gerono](https://en.wikipedia.org/wiki/Lemniscate_of_Gerono) and trains an SOM
+ * that finds that circular pattern. \param[out] data matrix to store data in
+ * \param[in] N number of points required
+ */
+void test_3d_classes(double *const *data, int N)
+{
+    const double R = 0.1; // radius of cluster
+    int i;
+    const int num_classes = 4;
+    const double centres[][3] = {
+        // centres of each class cluster
+        {.5, .5, .5},   // centre of class 1
+        {.5, -.5, -.5}, // centre of class 2
+        {-.5, .5, .5},  // centre of class 3
+        {-.5, -.5 - .5} // centre of class 4
+    };
+
+#ifdef _OPENMP
+#pragma omp for
+#endif
+    for (i = 0; i < N; i++)
+    {
+        int class = rand() % num_classes; // select a random class for the point
+
+        // create random coordinates (x,y,z) around the centre of the class
+        data[i][0] = _random(centres[class][0] - R, centres[class][0] + R);
+        data[i][1] = _random(centres[class][1] - R, centres[class][1] + R);
+        data[i][2] = _random(centres[class][2] - R, centres[class][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);
+        */
+    }
+}
+
+/** Test that creates a random set of points distributed in six clusters in
+ * 3D space. The following
+ * [CSV](https://en.wikipedia.org/wiki/Comma-separated_values) files are created
+ * to validate the execution:
+ * * `test3.csv`: random test samples points with a circular pattern
+ * * `w31.csv`: initial random map
+ * * `w32.csv`: trained SOM map
+ *
+ * The outputs can be readily plotted in [gnuplot](https:://gnuplot.info) using
+ * the following snippet
+ * ```gnuplot
+ * set datafile separator ','
+ * plot "test3.csv" title "original", \
+ *      "w31.csv" title "w1", \
+ *      "w32.csv" title "w2"
+ * ```
+ */
+void test3()
+{
+    int j, N = 200;
+    int features = 3;
+    int num_out = 20;
+    double **X = (double **)malloc(N * sizeof(double *));
+    double **W = (double **)malloc(num_out * sizeof(double *));
+    for (int i = 0; i < (num_out > N ? num_out : N); i++)
+    {
+        if (i < N) // only add new arrays if i < N
+            X[i] = (double *)malloc(features * sizeof(double));
+        if (i < num_out) // only add new arrays if i < num_out
+        {
+            W[i] = (double *)malloc(features * sizeof(double));
+
+#ifdef _OPENMP
+#pragma omp for
+#endif
+            // preallocate with random initial weights
+            for (j = 0; j < features; j++)
+                W[i][j] = _random(-1, 1);
+        }
+    }
+
+    test_3d_classes(X, N); // create test data around the lamniscate
+    save_nd_data("test3.csv", X, N, features); // save test data points
+    save_nd_data("w31.csv", W, num_out,
+                 features); // save initial random weights
+    kohonen_som_tracer(X, W, N, features, num_out, 0.01); // train the SOM
+    save_nd_data("w32.csv", W, num_out, features); // save the resultant weights
+
+    for (int i = 0; i < (num_out > N ? num_out : N); i++)
+    {
+        if (i < N)
+            free(X[i]);
+        if (i < num_out)
+            free(W[i]);
+    }
+    free(X);
+    free(W);
+}
+
+/**
+ * Convert clock cycle difference to time in seconds
+ *
+ * \param[in] start_t start clock
+ * \param[in] start_t end clock
+ * \returns time difference in seconds
+ */
+double get_clock_diff(clock_t start_t, clock_t end_t)
+{
+    return (double)(end_t - start_t) / (double)CLOCKS_PER_SEC;
+}
+
+/** Main function */
+int main(int argc, char **argv)
+{
+#ifdef _OPENMP
+    printf("Using OpenMP based parallelization\n");
+#else
+    printf("NOT using OpenMP based parallelization\n");
+#endif
+    clock_t start_clk = clock();
+    test1();
+    clock_t end_clk = clock();
+    printf("Test 1 completed in %.4g sec\n",
+           get_clock_diff(start_clk, end_clk));
+    start_clk = clock();
+    test2();
+    end_clk = clock();
+    printf("Test 2 completed in %.4g sec\n",
+           get_clock_diff(start_clk, end_clk));
+    start_clk = clock();
+    test3();
+    end_clk = clock();
+    printf("Test 3 completed in %.4g sec\n",
+           get_clock_diff(start_clk, end_clk));
+    printf("(Note: Calculated times include: creating test sets, training "
+           "model and writing files to disk.)\n\n");
+    return 0;
+}