Kohonen SOM trace/chain algorithm

Collaboration diagram for Kohonen SOM trace/chain algorithm:

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	_random (double a, double b)
	Helper function to generate a random number in a given interval.
int	save_nd_data (const char *fname, double **X, int num_points, int num_features)
	Save a given n-dimensional data martix to file.
void	kohonen_get_min_1d (double const *X, int N, double *val, int *idx)
	Get minimum value and index of the value in a vector.
void	kohonen_update_weights (double const *x, double *const *W, double *D, int num_out, int num_features, double alpha, int R)
	Update weights of the SOM using Kohonen algorithm.
void	kohonen_som_tracer (double **X, double *const *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

a	lower limit
b	upper limit

Returns

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

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

kohonen_get_min_1d()

void kohonen_get_min_1d	(	double const *	X,
		int	N,
		double *	val,
		int *	idx
	)

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

Parameters

[in]	X	vector to search
[in]	N	number of points in the vector
[out]	val	minimum value found
[out]	idx	index where minimum value was found

{
    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];
        }
    }
}

kohonen_som_tracer()

void kohonen_som_tracer	(	double **	X,
		double *const *	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 alpha = 1.f;
    double *D = (double *)malloc(num_out * sizeof(double));
 
    // Loop alpha from 1 to alpha_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
            kohonen_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);
}

Here is the call graph for this function:

kohonen_update_weights()

void kohonen_update_weights	(	double const *	x,
		double *const *	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

{
    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 smallest Euclidian distance to
    // the current pattern
    int d_min_idx;
    double d_min;
    kohonen_get_min_1d(D, num_out, &d_min, &d_min_idx);
 
    // step 3a: get the neighborhood range
    int from_node = max(0, d_min_idx - R);
    int to_node = 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++)
        for (k = 0; k < num_features; k++)
            // update weights of nodes in the neighborhood
            W[j][k] += alpha * (x[k] - W[j][k]);
}

Here is the call graph for this function:

save_nd_data()

int save_nd_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 overwriten 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
                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;
}