/* koho - Self-Organizing Feature Maps - (Kohonen Networks).
 *
 * Written by Jim Warhol, for CSC 671 Winter '95.
 *
 * koho uses Teuvo Kohonen's unsupervised learning neural network model
 * to "learn" the topology of a set of input patterns.
 *
 * Compile with: gcc koho.c -o bin/koho -lm
 *
 * execute with: koho loop_count
 *
 * loop_count is the number of times to iterate through all input patterns.
 *
 * The patterns file is assumed to be in the following format:
 * 3 50 0.15		(# of ideal points, # of input patterns, fuzz factor)
 * A  0.0  0.0  1.0	(First ideal data point, for results analysis)
 * B  1.0  0.0  0.0	(Second ideal data point)
 * C  0.0  0.0 -1.0	(Third ideal data point, followed by others if more)
 * x.xxx x.xxx x.xxx (based on x x.xxx x.xxx x.xxx)  (First pattern)
 * (repeated to equal as many input patterns as indicated)
 */

#include <math.h>
#include <stdio.h>

/*
 * global definitions.
 */

#define INITALPHA  0.30	/* initial winner adjustment factor */
#define INITBETA   0.10	/* initial neighbor adjustment factor */
#define LEARNLOSS  0.9	/* periodic reduction in learning adjustment */
#define IDEALMAX   10	/* maximum number of ideal data points */
#define PATTERNMAX 200	/* maximum number of input patterns */

int i, j;			/* scratch variables */

float alpha = INITALPHA;	/* winner learning adjustment factor */
float beta = INITBETA;		/* neighbor learning adjustment factor */
int loop_count;			/* number of desired iterations thru patterns */
int ideal_count;		/* number of ideal data points */
int pattern_count;		/* number of patterns in patterns file */
float fuzz_factor;		/* fuzz factor used in pattern generation */
int current_pattern;		/* current pattern being inspected */
int winner_x, winner_y;		/* grid coordinates of winning unit */

struct vector { float x; float y; float z; };

struct vector Ideal[IDEALMAX];
struct vector Weights[5][5];
struct vector Patterns[PATTERNMAX];
struct vector V;			/* scratch vector */

char Ideal_name[IDEALMAX][20];
char patterns_name[14] = "patterns";
FILE *patterns;

/*
 * function prototypes.
 */

void adjust_weights();
void adjust_unit_weights(int x, int y, float adjustment);
void analyze_results();
void calculate_greatest_activation();
int check_range(int x, int y);
float dot_product(struct vector V1, struct vector V2);
void initialize_koho(int argc, char *argv[]);
void initialize_weights();
void load_patterns();
struct vector normalize(struct vector V);
float pickrand();

/* main program */

void main(int argc, char *argv[])
{

   int cycle;

   initialize_koho(argc, argv); /* initialize patterns & random numbers */

   initialize_weights(); /* initialize Weights with random unit vectors */

   load_patterns();	 /* read input pattern data points from patterns */

   /* Cycle through all patterns loop_count times, training the network */

   for (cycle = 1; cycle <= loop_count; cycle++)
      for (current_pattern=0;current_pattern<pattern_count;current_pattern++) {
         calculate_greatest_activation();
         adjust_weights();

         /* every other cycle, reduce learning rates */

         if (cycle % 10 == 0) alpha *= LEARNLOSS, beta *= LEARNLOSS;
      }

   analyze_results();

   exit(0);

}

/* void adjust_unit_weights(int x, int y, float adjustment);
 * 
 * adjust_unit_weights adjusts the weights of grid unit (x,y) by adding
 * the current pattern coordinates times the specified adjustment.
 * No action is taken if either of the grid coordinates x or y are
 * out of range.
 */

void adjust_unit_weights(int x, int y, float adjustment)
{

   if (check_range(x, y)) {
      V.x = Weights[x][y].x + adjustment*Patterns[current_pattern].x;
      V.y = Weights[x][y].y + adjustment*Patterns[current_pattern].y;
      V.z = Weights[x][y].z + adjustment*Patterns[current_pattern].z;
      Weights[x][y] = normalize(V);
   }
   return;
}

/* void adjust_weights();
 *
 * adjust weights for winning unit by alpha and neighbors of winning
 * unit by beta. alpha and beta are learning rates, which decrease
 * over time.
 */

void adjust_weights()
{

   /* adjust weights of winning unit and neighboring units */

   adjust_unit_weights(winner_x - 1, winner_y + 1, beta);
   adjust_unit_weights(winner_x    , winner_y + 1, beta);
   adjust_unit_weights(winner_x + 1, winner_y + 1, beta); 

   adjust_unit_weights(winner_x - 1, winner_y    , beta); 
   adjust_unit_weights(winner_x    , winner_y    , alpha);
   adjust_unit_weights(winner_x + 1, winner_y    , beta); 

   adjust_unit_weights(winner_x - 1, winner_y - 1, beta); 
   adjust_unit_weights(winner_x    , winner_y - 1, beta); 
   adjust_unit_weights(winner_x + 1, winner_y - 1, beta); 

   return;
}

/* void analyze_results();
 *
 * Analyze results of Kohonen Network.
 */

void analyze_results()
{
   int loop;			/* loop counter for dot product comparison */
   int closest_ideal;		/* index of ideal point closest to weight */
   float dot;			/* dot product of ideal point with weight */
   float current_max;		/* current dot product closest to ideal point */

   printf("\nResults after %d iterations through %d patterns, totalling %d"
      " patterns.\n",
      loop_count, pattern_count, loop_count * pattern_count);
   printf("(Patterns were generated with a fuzz factor of %f.)\n", fuzz_factor);

   for (i = 0; i < 5; i++) {
      for (j = 0; j < 5; j++) {

         /* determine closest ideal data point to the unit */

         current_max = -999.9;
         for (loop = 0; loop < ideal_count; loop++)
            if ((dot = dot_product(Weights[i][j],Ideal[loop])) > current_max) {
               current_max = dot, closest_ideal = loop;
            }

         /* print out name of closest ideal data point for this unit */

         printf("%s ", Ideal_name[closest_ideal]);
      }
      printf("\n");
   }

   return;
}

/* void calculate_greatest_activation();
 */

void calculate_greatest_activation()
{

   float winning_activation;	/* activation sum for winner, so far */
   float this_activation;	/* activation sum for current grid unit */

   winner_x = 0; winner_y = 0;	/* initialize winner coordinates */

   winning_activation = -999.9;	/* initialize activation sum of winner */

   /* loop through the grid to find the unit with the greatest activation sum */

   for (i = 0; i < 5; i++)
      for (j = 0; j < 5; j++)
         if ((this_activation = dot_product(Patterns[current_pattern], 
            Weights[i][j])) > winning_activation)
               winning_activation = this_activation, winner_x = i, winner_y = j;

   return;
}

/* int check_range(int x, int y);
 */

int check_range(int x, int y)
{
   return x >= 0 && x <= 4 && y >= 0 && y <= 4;
}

/* float dot_product(struct vector V1, struct vector V2);
 */

float dot_product(struct vector V1, struct vector V2)
{
   return (V1.x * V2.x) + (V1.y * V2.y) + (V1.z * V2.z);
}

/* void initialize_koho(int argc, char *argv[]);
 * 
 * initialize_koho initializes the patterns patterns file and initializes
 * the random number generator.
 */

void initialize_koho(int argc, char *argv[])
{
   if (argc != 2) {
      fprintf(stderr, "ERROR: You must specify the loop count. e.g., "
         "koho 10\n");
      exit(1);
   }

   if (sscanf(argv[1], "%d", &loop_count) != 1) {
      fprintf(stderr, "ERROR: Loop count is in error.\n");
      exit(1);
   }
   
   /* open the patterns patterns file */

   if ((patterns = fopen(patterns_name, "r")) == NULL) {
      fprintf(stderr, "ERROR: cannot open %s.\n", patterns_name);
      exit(1); }

   /* read the number of ideal data points and the number of patterns */

   if (fscanf(patterns, "%d %d %f ", &ideal_count, &pattern_count,
      &fuzz_factor) != 3) {
      fprintf(stderr, "ERROR: fscanf returned other than 3.\n");
      exit(1); }

   if (ideal_count < 1 || ideal_count > IDEALMAX) {
      fprintf(stderr, "ERROR: Ideal point count in patterns is out "
         "of range (1 to %d).\n", IDEALMAX);
      exit(1); }

   if (pattern_count < 1 || pattern_count > PATTERNMAX) {
      fprintf(stderr, "ERROR: More than %d patterns defined in patterns.\n", 
         PATTERNMAX);
      exit(1); }

   /* read the name and value of the ideal data points */

   for (i = 0; i < ideal_count; i++)	/* read in ideal points */
      if (fscanf(patterns, "%s %f %f %f ",
         &Ideal_name[i], &Ideal[i].x, &Ideal[i].y, &Ideal[i].z) != 4) {
         fprintf(stderr, "ERROR: fscanf returned other than 4.\n");
         exit(1); }

   /* initialize random number generator */

   srandom((int) time()); 

}

/* void initialize_weights();
 */

void initialize_weights()
{

   for (i = 0; i < 5; i++)
      for (j = 0; j < 5; j++) {
         V.x = pickrand(); V.y = pickrand(); V.z = pickrand();
         Weights[i][j] = normalize(V);
      }
   return;

}

/* void load_patterns();
 *
 * load_patterns loads the previously generated patterns used to train
 * the Kohonen network from file "patterns". These patterns are loaded
 * into the Patterns array. Patterns are placed in the array in random
 * order.
 */

void load_patterns()
{
   int fscode;
   float x, y, z;

   /* initialize all pattern x coordinates with a flag indicating
      the pattern has not yet been loaded into that cell*/

   for (i = 0; i < pattern_count; i++) Patterns[i].x = -999.9;

   /* read in each pattern, and find a random (unused) spot in the 
      pattern array */

   for (i = 0; i < pattern_count; i++) {
      if (fscode = fscanf(patterns, "%f %f %f %*43c ", &x, &y, &z) != 3) {
         fprintf(stderr, "ERROR: fscanf returned %i; 3 expected.\n", fscode);
         exit(1); }

      /* find a suitable spot in the patterns array for this pattern */

      while (Patterns[j = (random() % pattern_count)].x > -999.0); 
      Patterns[j].x = x, Patterns[j].y = y, Patterns[j].z = z;
      }

   return;
}

/* struct vector normalize(struct vector V);
 *
 * Vnew = normalize(Vold);
 *
 * normalizes vector Vold so it has length 1 (and thus falls on the unit 
 * sphere).
 */

struct vector normalize(struct vector V)
{
   struct vector Result;
   double vector_length;
 
   vector_length = sqrt(V.x * V.x + V.y * V.y + V.z * V.z);

   Result.x = V.x / vector_length;
   Result.y = V.y / vector_length;
   Result.z = V.z / vector_length;

   return Result;
};

/* float pickrand()
 *
 * pickrand generates a random number in the range -1.0 to +1.0.
 */

float pickrand()
{
   return fmod(random(), 2.0 * 10000001.0) / 10000000.0 - 1.0;
}