common.c

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

  }
  if (x == sdim && sdim != data->ldim){
    /* No mapping info available (all struct info is (-1 -1)).
       Will encourage the generation of random mapping info... */
    for (x = data->ldim; x < sdim; x+=2){
      maval[x]   = map->xdim -1;
      maval[x+1] = map->ydim -1;
      mival[x]   = 0;
      mival[x+1] = 0;
    }
  }

  switch (mode){
  case INIT_LINEAR:
    fprintf(stderr, "\nlinear mode\n");
    /* Linear initialization of vector values */
    ComputeDistance = ComputeHexaDistance;

    nc = 0;
    maxdist = ComputeDistance(map->xdim-1, map->ydim-1, 0, 0);
    for (y = 0u; y < map->ydim; y++){
      for (x = 0u; x < map->xdim; x++){
	dist = ComputeDistance(map->codes[nc].x, map->codes[nc].y, 0, 0);
	dist = sqrtf((FLOAT)dist);
	for (i = 0; i < dim; i++){
	  map->codes[nc].points[i] = mival[i] + ((maval[i] - mival[i]) * dist / maxdist);
	}
	for (; i < map->dim; i++){ /* In VQ mode, map->dim > dim */
	  if (i-dim == nc)
	    map->codes[nc].points[i] = 1.0;
	  else
	    map->codes[nc].points[i] = 0.0;
	}
	nc++;
      }
    }
    break;
  case INIT_RANDOM:
  default:
    //    fprintf(stderr, "\nrandom mode\n");
    /* Randomize the vector values */
    for (i = 0; i < noc; i++){
      for (x = 0; x < dim; x++){
	map->codes[i].points[x] = mival[x] + (maval[x] - mival[x]) * drand48();
      }
      for (; x < map->dim; x++){ /* In VQ mode, map->dim > dim */
	if (x-dim == i)
	  map->codes[i].points[x] = 1.0;
	else
	  map->codes[i].points[x] = 0.0;
      }
    }
    break;
  }

  free(mival);
  free(maval);

  return 1;
}

/******************************************************************************
Description: This function is called in VQ mode, and will compute the auxillary
             variables a and b.

Return value: This function does not return a value.
******************************************************************************/
void VQSet_ab(struct Parameters *parameters)
{
  UNSIGNED n, i, noc;
  struct Map *map;
  struct Graph *graph;
  FLOAT tmp;

  map = &parameters->map;
  if ((graph = parameters->train) == NULL)
    if ((graph = parameters->valid) == NULL)
      if ((graph = parameters->test) == NULL){
	AddError("Init error in VQSet_ab: No data available.");
	return;
      }

  noc = map->ydim * map->xdim;
  for (n = 0; n < noc; n++){
    tmp = 0.0;
    for (i = 0; i < graph->FanOut * noc; i++)
      tmp += SQR(map->codes[n].points[graph->ldim + i]);
    map->codes[n].a = tmp;
    tmp = 0.0;
    for (i = 0; i < graph->FanIn * noc; i++)
      tmp += SQR(map->codes[n].points[graph->FanOut * noc + graph->ldim + i]);
    map->codes[n].b = tmp;
  }
}

/*****************************************************************************
Description: Compute optimal mu-weights

Return value: This function does not return a value.
*****************************************************************************/
void GetMuValues(struct Parameters *params, FLOAT *mu1, FLOAT *mu2, FLOAT *mu3, FLOAT *mu4)
{
  UNSIGNED ldim, cend, tend, pend;
  UNSIGNED dim, num;
  UNSIGNED i, x;
  FLOAT *avg, *sigma;
  FLOAT d1, d2, d3, d4;
  FLOAT x1, x2, x3, x4, k;
  struct Graph *gptr;
  struct Node *node;

  if (params->map.codes == NULL || params->map.dim == 0 || params->train==NULL)
    return;

  if(params->map.topology == TOPOL_VQ){
    fprintf(stderr, "\nWarning: Suggested mu values in VQ mode are incorrect.\n");
    fprintf(stderr, "         Function GetMuValues() not yet adapted to VQ mode.\n");
  }

  dim = params->map.dim;
  avg    = (FLOAT*)MyCalloc(dim, sizeof(FLOAT));
  sigma  = (FLOAT*)MyCalloc(dim, sizeof(FLOAT));

  /* Compute average values */
  ldim = params->train->ldim;
  cend = ldim + 2*params->train->FanOut;
  pend = cend + 2*params->train->FanIn;
  tend = pend + params->train->tdim;

  if(dim != tend){
    AddError("Dimension of codebooks != dimension of train-set vectors");
    *mu1 = 0;
    *mu2 = 0;
    *mu3 = 0;
    *mu4 = 0;
    return;
  }

  num = 0;
  for (gptr = params->train; gptr != NULL; gptr = gptr->next){
    for (i = 0; i < gptr->numnodes; i++){
      node = gptr->nodes[i];
      for (x = 0; x < ldim; x++)    /* Sum all values in node label   */
	avg[x] += node->points[x];
      for (x = pend; x < tend; x++) /* Sum all values in target vector */
	avg[x] += node->points[x];
      num++;
    }
  }
  for (i = 0; i < dim; i++)        /* Compute the average */
    avg[i] = avg[i]/num;
  for (x = ldim; x < pend; x+=2){  /* Average of parent and child states */
    avg[x]   = (params->map.xdim-1)/2.0;
    avg[x+1] = (params->map.ydim-1)/2.0;
  }

  /* Compute std.deviation */
  num = 0;
  for (gptr = params->train; gptr != NULL; gptr = gptr->next){
    for (i = 0; i < gptr->numnodes; i++){
      node = gptr->nodes[i];
      for (x = 0; x < ldim; x++)
	sigma[x] += SQR(node->points[x] - avg[x]);
      for (x = pend; x < tend; x++)
	sigma[x] += SQR(node->points[x] - avg[x]);
      num++;
    }
  }
  for (i = 0; i < dim; i++)
    sigma[i] = sigma[i]/num;
  for (x = ldim; x < pend; x++)
    sigma[x] = SQR(avg[x]/2);

  for (d1 = 0.0, x = 0; x < ldim; x++)
    d1 += sigma[x];
  for (d2 = 0.0; x < cend; x++)
    d2 += sigma[x];
  for (d3 = 0.0; x < pend; x++)
    d3 += sigma[x];
  for (d4 = 0.0; x < tend; x++)
    d4 += sigma[x];

  if (ldim > 0){              /* If we have a node label */
    x1 = 1.0;
    x2 = (d2 > 0.0) ? d1/d2 : 0.0;
    x3 = (d3 > 0.0) ? d1/d3 : 0.0;
    x4 = (d4 > 0.0) ? d1/d4 : 0.0;
    k = 1/(x1+x2+x3+x4); /* Normalizing factor */
  }
  else if (cend - ldim > 0){  /* No node label but child states */
    x1 = 0.0;
    x2 = 1.0;
    x3 = (d3 > 0.0) ? d2/d3 : 0.0;
    x4 = (d4 > 0.0) ? d2/d4 : 0.0;
    k = 1/(x2+x3+x4);    /* Normalizing factor */
  }
  else if (pend - cend > 0){  /* No label or child states but parent states */
    x1 = 0.0;
    x2 = 0.0;
    x3 = 1.0;
    x4 = (d4 > 0.0) ? d3/d4 : 0.0;
    k = 1/(x3+x4);       /* Normalizing factor */
  }
  else if (tend - pend > 0){  /* Only have target vector */
    x1 = 0.0;
    x2 = 0.0;
    x3 = 0.0;
    x4 = 1.0;
    k = 1.0;
  }
  else{                       /* Node dimension is empty */
    x1 = 0.0;
    x2 = 0.0;
    x3 = 0.0;
    x4 = 0.0;
    k = 0.0;
  }

  *mu1 = k*x1;
  *mu2 = k*x2;
  *mu3 = k*x3;
  *mu4 = k*x4;

  free(avg);
  free(sigma);
}

/*****************************************************************************
Description: Suggest optimal mu-weights. Suggestions are printed to stderr.

Return value: This function does not return a value.
*****************************************************************************/
void SuggestMu(struct Parameters *params)
{
  FLOAT mu1, mu2, mu3, mu4;
  char buffer[128];

  GetMuValues(params, &mu1, &mu2, &mu3, &mu4);

  if (similar(params->mu1 + params->mu2 + params->mu3 + params->mu4, 0.0, 4.0*EPSILON)){
    params->mu1 = mu1;
    params->mu2 = mu2;
    params->mu3 = mu3;
    params->mu4 = mu4;
    strcpy(buffer, "Will use:");
    if (mu1 > 0.00001)
      sprintf(&buffer[strlen(buffer)], " -mu1 %.9f", mu1);
    else if (mu1 > 0.0)
      sprintf(&buffer[strlen(buffer)], " -mu1 %.9e", mu1);
    if (mu2 > 0.00001)
	sprintf(&buffer[strlen(buffer)], " -mu2 %.9f", mu2);
    else if (mu4 > 0.0)
	sprintf(&buffer[strlen(buffer)], " -mu2 %.9e", mu2);
    if (mu3 > 0.00001)
	sprintf(&buffer[strlen(buffer)], " -mu3 %.9f", mu3);
    else if (mu4 > 0.0)
	sprintf(&buffer[strlen(buffer)], " -mu3 %.9e", mu3);
    if (mu4 > 0.00001)
	sprintf(&buffer[strlen(buffer)], " -mu4 %.9f", mu4);
    else if (mu4 > 0.0)
	sprintf(&buffer[strlen(buffer)], " -mu4 %.9e", mu4);
    AddMessage(buffer);
  }
  else{
    if (!approx(params->mu1, mu1, 0.01) || !approx(params->mu2, mu2, 0.01) || !approx(params->mu3, mu3, 0.01) || !approx(params->mu4, mu4, 0.01)){
      AddMessage("Caution: The mu-values are not optimal.");
      strcpy(buffer, "Suggesting use of:");
      if (mu1 > 0.00001)
	sprintf(&buffer[strlen(buffer)], " -mu1 %1.9f", mu1);
      else if (mu1 > 0.0)
	sprintf(&buffer[strlen(buffer)], " -mu1 %.9e", mu1);
      if (mu2 > 0.00001)
	sprintf(&buffer[strlen(buffer)], " -mu2 %1.9f", mu2);
      else if (mu2 > 0.0)
	sprintf(&buffer[strlen(buffer)], " -mu2 %.9e", mu2);
      if (mu3 > 0.00001)
	sprintf(&buffer[strlen(buffer)], " -mu3 %1.9f", mu3);
      else if (mu3 > 0.0)
	sprintf(&buffer[strlen(buffer)], " -mu3 %.9e", mu3);
      if (mu4 > 0.00001)
	sprintf(&buffer[strlen(buffer)], " -mu4 %1.9f", mu4);
      else if (mu4 > 0.0)
	sprintf(&buffer[strlen(buffer)], " -mu4 %.9e", mu4);
      AddMessage(buffer);
    }
  }
}