testsom.c

SOM-sd 是现在国外非常流行的一个神经网络的结构自组织特征映射算法

/*
  THIS IS AN UNTESTED DEVELOPMENT VERSION!
  -> NOT FOR DISTRIBUTION!
  -> USE IT WITH CARE BUT NOTE THAT
     if testsom doesn't meet your expectations then write your own code!

  Contents: The main module used to evaluate a somsd

  Author: Markus Hagenbuchner
          [e-mail]
 */


/************/
/* Includes */
/************/
#include <ctype.h>
#include <limits.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include "common.h"
#include "data.h"
#include "fileio.h"
#include "system.h"
#include "train.h"
#include "utils.h"

/*************/
/* Constants */
/*************/
#define CONTEXT        0x00000001
#define CONTEXTUAL     0x00000002
#define MAPNODETYPE    0x00000004
#define MAPNODELABEL   0x00000008
#define MAPGRAPH       0x00000010
#define MAPNODE        0x00000020
#define SHOWGRAPH      0x00000040
#define SHOWCLUSTERING 0x00000080
#define SHOWSUBGRAPHS  0x00000100
#define PRECISION      0x00000200
#define RETRIEVALPERF  0x00000400
#define CLASSIFY       0x00000800
#define ANALYSE        0x00001000
#define BALANCE        0x00002000
#define DISTANCES      0x00004000
#define WEBSOM         0x00008000
#define TRUNCATE       0x00010000
#define QUIET          0x01000000

struct VMap{
  UNSIGNED max;
  UNSIGNED **activation;
  UNSIGNED numclasses;
  UNSIGNED ***classes;
  UNSIGNED **winnerclass;
};

struct AllHits{
  struct Graph *graph;
  struct Node *node;
  char *structID;
  char *substructID;
  struct AllHits *next;
};

int KstepEnabled = 0;

/* Begin functions... */

/******************************************************************************
Description: Print usage information to the screen

Return value: The function does not return.
******************************************************************************/
void Usage()
{
  fprintf(stderr, "\n\
Usage: testsom [options]\n\n\
Options are:\n\
    -cin <fname>        Codebook file\n\
    -din <fname>        The file which holds the training data set.\n\
    -tin <fname>        The file which holds the test data set.\n\
    -mode <mode>        Test mode, which can be:\n\
                        context   Writes a data set suitable for further\n\
                                  contextual processing (multi-layer contexttual mode).\n\
                        contextual Computes stable point of dataset, then writes\n\
                                  the result to stdout (single layer contextual mode).\n\
                        mapgraph  For each codebook, plot (at most) one graph\n\
                                  structure that was mapped at that location.\n\
                        mapnode   For each codebook, list all nodes that were\n\
                                  mapped at that location.\n\
                        mapcluster Same as mapgraph, but plot class membership as well.\n\
                        showgraph Visualize a graph at a given location\n\
                                  location info is read from stdin.\n\
                        showsubgraph Visualize a given graph and all its\n\
                                  subgraphs. info is read from stdin.\n\
                        maplabel  For each codebook, list the labels of root\n\
                                  nodes that activated the codebook.\n\
                        precision Compute the mapping precision of the map.\n\
                        retrievalperformance Compute the retrieval\n\
                                  performance.\n\
                        classify  Classify a given set of data. A labelled\n\
                                  training set needs to be available.\n\
                        analyse   Statistically analyse a given dataset.\n\
                        balance   Produce a balanced dataset\n\
                        truncate <n>  Truncate outdegree of graphs in a given\n\
                                  dataset to at most n.\n\
                        distances List for every node in the dataset distances\n\
                                  to every codebook in the map.\n\
    -mu1 float[:float]  Weight(range) for the label component.\n\
    -mu2 float[:float]  Weight(range) for the child state component.\n\
    -mu3 float[:float]  Weight(range) for the parents position component.\n\
    -mu4 float[:float]  Weight(range) for the class label component.\n\
    -quiet              Restrict amount of text printed to screen.\n\
    -help               Print this help.\n\
 \n");
  exit(0);
}

/*****************************************************************************
Description: Print an xfig header

Return value:
*****************************************************************************/
void PrintXfigHeader(FILE *ofile)
{
  time_t t;

  t = clock();
  fprintf(ofile, "#FIG 3.2\n");
  fprintf(ofile, "Landscape\n");
  fprintf(ofile, "Center\n");
  fprintf(ofile, "Inches\n");
  fprintf(ofile, "Letter\n");
  fprintf(ofile, "100.00\n");
  fprintf(ofile, "Single\n");
  fprintf(ofile, "-2\n");
  fprintf(ofile, "1200 2\n");
  fprintf(ofile, "0 32 #bebabe\n");

  fprintf(ofile, "0 33 #f6f6f6\n");
  fprintf(ofile, "0 34 #f0f0f0\n");
  fprintf(ofile, "0 35 #e6e6e6\n");
  fprintf(ofile, "0 36 #e0e0e0\n");
  fprintf(ofile, "0 37 #d6d6d6\n");
  fprintf(ofile, "0 38 #d0d0d0\n");
  fprintf(ofile, "0 39 #c6c6c6\n");
  fprintf(ofile, "0 40 #c0c0c0\n");
  fprintf(ofile, "0 41 #b6b6b6\n");
  fprintf(ofile, "0 42 #b0b0b0\n");
  fprintf(ofile, "0 43 #a6a6a6\n");
  fprintf(ofile, "0 44 #a0a0a0\n");
  fprintf(ofile, "0 45 #969696\n");
  fprintf(ofile, "0 46 #909090\n");
}

void DrawPattern(int pattern, int hoff, int voff, int w, int h)
{
  pattern = 15;

  switch (pattern){
  case 0:
    printf("2 1 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 2\n");
    printf("\t%d %d %d %d\n", hoff+w/3, voff-h/2, hoff+2*w/3, voff+h/2);
    break;
  case 1:
    printf("2 1 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 2\n");
    printf("\t%d %d %d %d\n", hoff+2*w/3, voff-h/2, hoff+w/3, voff+h/2);
    break;
  case 2:
    printf("2 1 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 2\n");
    printf("\t%d %d %d %d\n", hoff+w/3, voff-h/2, hoff+2*w/3, voff+h/2);
    printf("2 1 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 2\n");
    printf("\t%d %d %d %d\n", hoff+2*w/3, voff-h/2, hoff+w/3, voff+h/2);
    break;
  case 3:
    printf("2 1 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 2\n");
    printf("\t%d %d %d %d\n", hoff+w/2, voff-h/2, hoff+w/2, voff+h/2);
    break;
  case 4:
    printf("2 1 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 2\n");
    printf("\t%d %d %d %d\n", hoff, voff, hoff+w, voff);
    break;
  case 5:
    printf("2 1 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 2\n");
    printf("\t%d %d %d %d\n", hoff+w/2, voff-h/2, hoff+w/2, voff+h/2);
    printf("2 1 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 2\n");
    printf("\t%d %d %d %d\n", hoff, voff, hoff+w, voff);
    break;
  case 6:
    printf("2 1 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 2\n");
    printf("\t%d %d %d %d\n", hoff+w/3, voff-h/2, hoff+2*w/3, voff+h/2);
    printf("2 1 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 2\n");
    printf("\t%d %d %d %d\n", hoff+2*w/3, voff-h/2, hoff+w/3, voff+h/2);
    printf("2 1 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 2\n");
    printf("\t%d %d %d %d\n", hoff+w/2, voff-h/2, hoff+w/2, voff+h/2);
    printf("2 1 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 2\n");
    printf("\t%d %d %d %d\n", hoff, voff, hoff+w, voff);
    break;
  case 7:
    printf("2 1 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 2\n");
    printf("\t%d %d %d %d\n", hoff+w/3, voff-h/2, hoff+w/3, voff+h/2);
    printf("2 1 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 2\n");
    printf("\t%d %d %d %d\n", hoff+w/2, voff-h/2, hoff+w/2, voff+h/2);
    printf("2 1 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 2\n");
    printf("\t%d %d %d %d\n", hoff+2*w/3, voff-h/2, hoff+2*w/3, voff+h/2);
    break;
  case 8:
    printf("2 1 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 4\n");
    printf("\t%d %d %d %d %d %d %d %d\n", hoff, voff, hoff+2*w/3, voff-h/2, hoff+2*w/3, voff+h/2, hoff, voff);
    break;
  case 9:
    printf("2 1 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 4\n");
    printf("\t%d %d %d %d %d %d %d %d\n", hoff+w/3, voff-h/2, hoff+w, voff, hoff+w/3, voff+h/2, hoff+w/3, voff-h/2);
    break;
  case 10:
    printf("2 1 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 2\n");
    printf("\t%d %d %d %d\n", hoff+w/3, voff-h/2, hoff+w/2, voff);
    printf("2 1 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 2\n");
    printf("\t%d %d %d %d\n", hoff+w/3, voff+h/2, hoff+w/2, voff);
    printf("2 1 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 2\n");
    printf("\t%d %d %d %d\n", hoff+w/2, voff, hoff+w, voff);
    break;
  case 11:
    printf("2 1 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 2\n");
    printf("\t%d %d %d %d\n", hoff, voff, hoff+w/2, voff);
    printf("2 1 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 2\n");
    printf("\t%d %d %d %d\n", hoff+w/2, voff, hoff+2*w/3, voff-h/2);
    printf("2 1 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 2\n");
    printf("\t%d %d %d %d\n", hoff+w/2, voff, hoff+2*w/3, voff+h/2);   
    break;
  case 12:
    printf("2 1 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 2\n");
    printf("\t%d %d %d %d\n", hoff+w/3, voff-h/2, hoff+w/2, voff);
    printf("2 1 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 2\n");
    printf("\t%d %d %d %d\n", hoff+2*w/3, voff-h/2, hoff+w/2, voff);
    printf("2 1 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 2\n");
    printf("\t%d %d %d %d\n", hoff+w/2, voff, hoff+w/2, voff+h/2);   
    break;
  case 13:
    printf("2 1 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 2\n");
    printf("\t%d %d %d %d\n", hoff+w/3, voff+h/2, hoff+w/2, voff);
    printf("2 1 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 2\n");
    printf("\t%d %d %d %d\n", hoff+2*w/3, voff+h/2, hoff+w/2, voff);
    printf("2 1 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 2\n");
    printf("\t%d %d %d %d\n", hoff+w/2, voff, hoff+w/2, voff-h/2);   
    break;
  case 14:
    printf("2 1 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 2\n");
    printf("\t%d %d %d %d\n", hoff, voff, hoff+w/2, voff);
    printf("2 1 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 2\n");
    printf("\t%d %d %d %d\n", hoff, voff, hoff+2*w/3, voff-h/2);
    printf("2 1 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 2\n");
    printf("\t%d %d %d %d\n", hoff, voff, hoff+2*w/3, voff+h/2);   
    break;
  case 15:
    printf("2 1 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 2\n");
    printf("\t%d %d %d %d\n", hoff+w/3, voff-h/2, hoff+w, voff);
    printf("2 1 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 2\n");
    printf("\t%d %d %d %d\n", hoff+w/3, voff+h/2, hoff+w, voff);
    printf("2 1 0 1 0 7 50 -1 -1 0.000 0 0 -1 0 0 2\n");
    printf("\t%d %d %d %d\n", hoff+w/2, voff, hoff+w, voff);
    break;
  }
}

/*****************************************************************************
Description: Draw the neurons of a map as .fig format

Return value:
*****************************************************************************/
void DrawMap(struct Map map, FLOAT scale, struct VMap *vmap)
{
  int w, h;       /* Width and height of a neuron */
  int hoff, voff; /* horizontal and vertical offsets */
  int x, y;       /* count through width and height of map */
  FLOAT cval;

  w = scale*1.3;
  h = scale*1.3;
  for (y = 0; y < map.ydim; y++){
    for (x = 0; x < map.xdim; x++){
      hoff = x*2*w/3 ;
      if (x % 2)
	voff = y*h-h/2;
      else
	voff = y*h;
      if (vmap != NULL){
	if (vmap->activation[y][x] == 0)
	  printf("2 1 0 2 32 7 50 -1 -1 0.000 0 0 -1 0 0 7\n");
	else{
	  cval = 46 - (vmap->max - vmap->activation[y][x])*13/vmap->max;
	  printf("2 1 0 2 32 %d 51 -1 20 0.000 0 0 -1 0 0 7\n", (int)cval);
	}
      }
      else
	printf("2 1 0 2 32 7 50 -1 -1 0.000 0 0 -1 0 0 7\n");

      printf("\t%d %d ", hoff, voff);
      printf("%d %d ", hoff+w/3, voff-h/2);
      printf("%d %d ", hoff+2*w/3, voff-h/2);
      printf("%d %d ", hoff+w, voff);
      printf("%d %d ", hoff+2*w/3, voff+h/2);
      printf("%d %d ", hoff+w/3, voff+h/2);
      printf("%d %d\n", hoff, voff);

      if (vmap->winnerclass && vmap->winnerclass[y][x] > 0)
	DrawPattern(vmap->winnerclass[y][x]-1, hoff, voff, w, h);
    }
  }
}

int GetGraphWidth(struct Graph* gptr)
{
  int n, max, maxdepth;
  int *stats;

  if (gptr == NULL)
    return 0;

  maxdepth = 0;
  for (n = 0; n < gptr->numnodes; n++){
    if (maxdepth < gptr->nodes[n]->depth)
      maxdepth = gptr->nodes[n]->depth;
  }
  stats = (int*)MyCalloc((maxdepth+1), sizeof(int));

  for (n = 0; n < gptr->numnodes; n++)
    stats[gptr->nodes[n]->depth]++;

  max = 0;
  for (n = 0; n < maxdepth+1; n++){
    if (max < stats[n])
      max = stats[n];
  }
  free(stats);
  return max;
}

void DrawXfigGraphAt(FILE *ofile, struct Graph* gptr, struct Node *root, int x, int y, float scale)
{
  struct Helper{
    int x, y;
    struct Node *node;
  };
//  static int offset = 11;
  int width, height;
  int d, cd, flag;
  int n, i, c, cn, j, wo, num;
  //  char buffer[1024];
  struct Helper **gmatrix, *xynode, *child;
  struct Node *node;
  //  int fontsize = 8;

  if (ofile == NULL || gptr == NULL)
    return;

  //Initializations
  width = GetGraphWidth(gptr);
  height = root->depth+1; //GetGraphDepth(gptr)+1;

  gmatrix = (struct Helper**)MyMalloc(height * sizeof(struct Helper*));
  for (i = 0; i < height; i++)
    gmatrix[i] = (struct Helper*)MyCalloc((width+1), sizeof(struct Helper));

  gmatrix[root->depth][0].node = root;
  for (d = root->depth; d > 0; d--){
    c = 0;
    for (i = 0; i < width; i++){
      if (gmatrix[d][i].node == NULL)
	continue;

      for (j = 0; j < gptr->FanOut; j++){
	if (gmatrix[d][i].node->children[j] != NULL){
	  gmatrix[d-1][c].node = gmatrix[d][i].node->children[j];
	  c++;
	}
      }
    }
  }
  for (d = 0; d < height; d++){
    num = 0;
    for (c = 0; c < width; c++){
      if (gmatrix[d][c].node != NULL)
	num++;
      gmatrix[d][c].x = scale+x+c*scale*3;
      gmatrix[d][c].y = y+(height-d)*scale*3;
    }
    wo = (3*(width-1-num)*scale)/2;
    for (c = 0; c < width; c++){
      gmatrix[d][c].x += wo;
    }
  }

  /*  for(n = 0; n < gptr->numnodes; n++){
    d = gptr->nodes[n]->depth;
    for (i = 0;gmatrix[d][i].node != NULL && gmatrix[d][i].node != gptr->nodes[n];i++);
    if (gmatrix[d][i].node == gptr->nodes[n])
      continue;
    gmatrix[d][i].node = gptr->nodes[n];
    gmatrix[d][i].x = i * spacing;
    gmatrix[d][i].y = 30 + d * spacing;
  }
  for (d = 0; d < height; d++){
    for(n = 0; gmatrix[d][n].node != NULL; n++);
    for(i = 0; gmatrix[d][i].node != NULL; i++){
      gmatrix[d][i].x = (spacing*(width-n+2*i))/2+offset;
    }
  }
  */

  flag = 0;
  for (d = 0; d < height; d++){
    for(n = 0; gmatrix[d][n].node != NULL; n++){
      xynode = &gmatrix[d][n];
      node = gmatrix[d][n].node;
      
      /* Print graph name */
      /*
      if (!flag && node->gname != NULL){
	fprintf(ofile, "4 0 0 50 0 0 %d 0.0000 4 105 645 ", (fontsize+4));
	fprintf(ofile, "%d %d %s\\001\n", offset*20, (fontsize+4)*10, node->gname);
	flag = 1;
      }
      */
      /* Print node target */
      /*
      if (node->lab.label > 1){
	sprintf(buffer, "c=\"%s\"", find_conv_to_lab(node->lab.label));
	fprintf(ofile, "4 0 0 50 0 0 %d 0.0000 4 105 645 ", fontsize);
	fprintf(ofile, "%d %d %s\\001\n", xynode->x*20 + 200+30, (int)(xynode->y*20-200+spacing+(fontsize*1.5)*10), buffer);
      }
      */
      /* Print node label */
      /*
      if (data->dimension > 0){
	if (node->points[0] == (int)node->points[0])
	  sprintf(buffer, "l=(%d", (int)node->points[0]);
	else
	  sprintf(buffer, "l=(%f", node->points[0]);
      }
      else
	strcpy(buffer, "(");
      */
      /*
      for (i = 1; i < data->dimension; i++){
	if (node->points[i] == (int)node->points[i])
	  sprintf(buffer, "%s, %d", buffer, (int)node->points[i]);
	else
	  sprintf(buffer, "%s, %f", buffer, node->points[i]);
      }
      strcat(buffer, ")");