Kohonen SOM topology algorithm

Collaboration diagram for Kohonen SOM topology algorithm:

Data Structures
struct	kohonen_array_3d
	to store info regarding 3D arrays More...

Macros
#define	max(a, b)   (((a) > (b)) ? (a) : (b))
	shorthand for maximum value
#define	min(a, b)   (((a) < (b)) ? (a) : (b))
	shorthand for minimum value

Functions
double *	kohonen_data_3d (const struct kohonen_array_3d *arr, int x, int y, int z)
	Function that returns the pointer to (x, y, z) ^th location in the linear 3D array given by:
double	_random (double a, double b)
	Helper function to generate a random number in a given interval.
int	save_2d_data (const char *fname, double **X, int num_points, int num_features)
	Save a given n-dimensional data martix to file.
int	save_u_matrix (const char *fname, struct kohonen_array_3d *W)
	Create the distance matrix or U-matrix from the trained weights and save to disk.
void	get_min_2d (double **X, int N, double *val, int *x_idx, int *y_idx)
	Get minimum value and index of the value in a matrix.
double	kohonen_update_weights (const double *X, struct kohonen_array_3d *W, double **D, int num_out, int num_features, double alpha, int R)
	Update weights of the SOM using Kohonen algorithm.
void	kohonen_som (double **X, struct kohonen_array_3d *W, int num_samples, int num_features, int num_out, double alpha_min)
	Apply incremental algorithm with updating neighborhood and learning rates on all samples in the given datset.

Detailed Description

Function Documentation

_random()

double _random	(	double	a,
		double	b
	)

Helper function to generate a random number in a given interval.

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 \([a,b)\)

   \[ y = (b - a) \times \frac{\text{(random number between 0 and RAND_MAX)} \; \text{mod}\; 100}{100} + a \]

Parameters

[in]	a	lower limit
[in]	b	upper limit

Returns

random number in the range \([a,b)\)

{
    return ((b - a) * (rand() % 100) / 100.f) + a;
}

get_min_2d()

void get_min_2d	(	double **	X,
		int	N,
		double *	val,
		int *	x_idx,
		int *	y_idx
	)

Get minimum value and index of the value in a matrix.

Parameters

[in]	X	matrix to search
[in]	N	number of points in the vector
[out]	val	minimum value found
[out]	x_idx	x-index where minimum value was found
[out]	y_idx	y-index where minimum value was found

{
    val[0] = INFINITY;  // initial min value
 
    for (int i = 0; i < N; i++)  // traverse each x-index
    {
        for (int j = 0; j < N; j++)  // traverse each y-index
        {
            if (X[i][j] < val[0])  // if a lower value is found
            {                      // save the value and its index
                x_idx[0] = i;
                y_idx[0] = j;
                val[0] = X[i][j];
            }
        }
    }
}

kohonen_data_3d()

double * kohonen_data_3d	(	const struct kohonen_array_3d *	arr,
		int	x,
		int	y,
		int	z
	)

Function that returns the pointer to (x, y, z) ^th location in the linear 3D array given by:

\[ X_{i,j,k} = i\times M\times N + j\times N + k \]

where \(L\), \(M\) and \(N\) are the 3D matrix dimensions.

Parameters

[in]	arr	pointer to kohonen_array_3d structure
[in]	x	first index
[in]	y	second index
[in]	z	third index

Returns

pointer to (x,y,z)^th location of data

{
    int offset = (x * arr->dim2 * arr->dim3) + (y * arr->dim3) + z;
    return arr->data + offset;
}

kohonen_som()

void kohonen_som	(	double **	X,
		struct kohonen_array_3d *	W,
		int	num_samples,
		int	num_features,
		int	num_out,
		double	alpha_min
	)

Apply incremental algorithm with updating neighborhood and learning rates on all samples in the given datset.

Parameters

[in]	X	data set
[in,out]	W	weights matrix
[in]	num_samples	number of output points
[in]	num_features	number of features per input sample
[in]	num_out	number of output points
[in]	alpha_min	terminal value of alpha

{
    int R = num_out >> 2, iter = 0;
    double **D = (double **)malloc(num_out * sizeof(double *));
    for (int i = 0; i < num_out; i++)
        D[i] = (double *)malloc(num_out * sizeof(double));
 
    double dmin = 1.f;  // average minimum distance of all samples
 
    // Loop alpha from 1 to slpha_min
    for (double alpha = 1.f; alpha > alpha_min && dmin > 1e-3;
         alpha -= 0.001, iter++)
    {
        dmin = 0.f;
        // 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 += kohonen_update_weights(X[sample], W, D, num_out,
                                           num_features, alpha, R);
        }
 
        // every 20th iteration, reduce the neighborhood range
        if (iter % 100 == 0 && R > 1)
            R--;
 
        dmin /= num_samples;
        printf("iter: %5d\t alpha: %.4g\t R: %d\td_min: %.4g\r", iter, alpha, R,
               dmin);
    }
    putchar('\n');
 
    for (int i = 0; i < num_out; i++) free(D[i]);
    free(D);
}

Here is the call graph for this function:

kohonen_update_weights()

double kohonen_update_weights	(	const double *	X,
		struct kohonen_array_3d *	W,
		double **	D,
		int	num_out,
		int	num_features,
		double	alpha,
		int	R
	)

Update weights of the SOM using Kohonen algorithm.

Parameters

[in]	X	data point
[in,out]	W	weights matrix
[in,out]	D	temporary vector to store distances
[in]	num_out	number of output points
[in]	num_features	number of features per input sample
[in]	alpha	learning rate \(0<\alpha\le1\)
[in]	R	neighborhood range

Returns

minimum distance of sample and trained weights

{
    int x, y, k;
    double d_min = 0.f;
 
#ifdef _OPENMP
#pragma omp for
#endif
    // step 1: for each 2D output point
    for (x = 0; x < num_out; x++)
    {
        for (y = 0; y < num_out; y++)
        {
            D[x][y] = 0.f;
            // compute Euclidian distance of each output
            // point from the current sample
            for (k = 0; k < num_features; k++)
            {
                double *w = kohonen_data_3d(W, x, y, k);
                D[x][y] += (w[0] - X[k]) * (w[0] - X[k]);
            }
            D[x][y] = sqrt(D[x][y]);
        }
    }
 
    // step 2:  get closest node i.e., node with smallest Euclidian distance to
    // the current pattern
    int d_min_x, d_min_y;
    get_min_2d(D, num_out, &d_min, &d_min_x, &d_min_y);
 
    // step 3a: get the neighborhood range
    int from_x = max(0, d_min_x - R);
    int to_x = min(num_out, d_min_x + R + 1);
    int from_y = max(0, d_min_y - R);
    int to_y = min(num_out, 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 = exp(-d2 / (2.f * alpha * alpha));
 
            for (k = 0; k < num_features; k++)
            {
                double *w = kohonen_data_3d(W, x, y, k);
                // update weights of nodes in the neighborhood
                w[0] += alpha * scale_factor * (X[k] - w[0]);
            }
        }
    }
    return d_min;
}

Here is the call graph for this function:

save_2d_data()

int save_2d_data	(	const char *	fname,
		double **	X,
		int	num_points,
		int	num_features
	)

Save a given n-dimensional data martix to file.

Parameters

[in]	fname	filename to save in (gets overwritten without confirmation)
[in]	X	matrix to save
[in]	num_points	rows in the matrix = number of points
[in]	num_features	columns in the matrix = dimensions of points

Returns

0 if all ok

-1 if file creation failed

{
    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
                fputc(',', fp);            // suffix comma
        }
        if (i < num_points - 1)  // if not the last row
            fputc('\n', fp);     // start a new line
    }
    fclose(fp);
    return 0;
}

save_u_matrix()

int save_u_matrix	(	const char *	fname,
		struct kohonen_array_3d *	W
	)

Create the distance matrix or U-matrix from the trained weights and save to disk.

Parameters

[in]	fname	filename to save in (gets overwriten without confirmation)
[in]	W	model matrix to save

Returns

0 if all ok

-1 if file creation failed

{
    FILE *fp = fopen(fname, "wt");
    if (!fp)  // error with fopen
    {
        char msg[120];
        sprintf(msg, "File error (%s): ", fname);
        perror(msg);
        return -1;
    }
 
    int R = max(W->dim1 >> 3, 2); /* neighborhood range */
 
    for (int i = 0; i < W->dim1; i++)  // for each x
    {
        for (int j = 0; j < W->dim2; j++)  // for each y
        {
            double distance = 0.f;
            int k;
 
            int from_x = max(0, i - R);
            int to_x = min(W->dim1, i + R + 1);
            int from_y = max(0, j - R);
            int to_y = min(W->dim2, j + R + 1);
            int l;
#ifdef _OPENMP
#pragma omp parallel for reduction(+ : distance)
#endif
            for (l = from_x; l < to_x; l++)  // scan neighborhoor in x
            {
                for (int m = from_y; m < to_y; m++)  // scan neighborhood in y
                {
                    double d = 0.f;
                    for (k = 0; k < W->dim3; k++)  // for each feature
                    {
                        double *w1 = kohonen_data_3d(W, i, j, k);
                        double *w2 = kohonen_data_3d(W, l, m, k);
                        d += (w1[0] - w2[0]) * (w1[0] - w2[0]);
                        // distance += w1[0] * w1[0];
                    }
                    distance += sqrt(d);
                    // distance += d;
                }
            }
 
            distance /= R * R;              // mean distance from neighbors
            fprintf(fp, "%.4g", distance);  // print the mean separation
            if (j < W->dim2 - 1)            // if not the last column
                fputc(',', fp);             // suffix comma
        }
        if (i < W->dim1 - 1)  // if not the last row
            fputc('\n', fp);  // start a new line
    }
    fclose(fp);
    return 0;
}

Here is the call graph for this function: