bayes.cc

KNN algorithm, and related

}
 
/* ----------------------------- MNI Header -----------------------------------
@NAME       : scale_matrices
@INPUT      : 
@OUTPUT     : 
@RETURNS    : 
@DESCRIPTION: scale matrices by given factor
@METHOD     : 
@GLOBALS    : num_classes, num_features
@CALLS      : 
@CREATED    : August 21, 1996 (John Sled)
@MODIFIED   : 
---------------------------------------------------------------------------- */
void scale_matrices(Real ***matrix, Real scale)
{
  int i, j, k;

  if(matrix) {
    for_less(i, 0, num_classes)
      for_less(j, 0, num_features)
        for_less(k, 0, num_features)
          matrix[i][j][k] *= scale;
  }
}


/* ----------------------------- MNI Header -----------------------------------
@NAME       : bayesian_classify_sample
@INPUT      : 
@OUTPUT     : 
@RETURNS    : sample class
@DESCRIPTION: given a feature vector and its size, return a class
@METHOD     : compute   Wi  where  Wi = Pi / sum(Pi)
                 Pi = (det(Si)/(2 pi)^d)^0.5 exp(-0.5 (x-u) Si (x -u)) p(i)
                 where
                    Si is the inverse of the ith covariance matrix
                    d  is the number of features
                    u  is the class mean vector
                    x  is the observed feature vector
                    p(i) is the apriori probability that tissue is ith class
@GLOBALS    : feature_vector, apriori
@CALLS      : 
@CREATED    : August 21, 1996 (John Sled)
@MODIFIED   : 
---------------------------------------------------------------------------- */
void bayesian_classify_sample(int *class_num, double *class_prob, 
                              int *class_labels)
{

  int           i, j, k;               /* counters */
  Real          exponent, sum, max_prob;
  int           class_index;

  /* calculate distance vectors */
  for_less ( i, 0, num_classes ) 
    for_less ( j, 0, num_features) 
      distance_vector[i][j] = feature_vector[j] - mean_feature_matrix[i][j];

  /* caluculate weights for each class */
  sum = 0.0;  /* sum of weights */
  for_less (i, 0, num_classes ) {
    
    exponent = 0.0;
    
    /* do diagonal */
    for_less(j, 0, num_features) {
      exponent += distance_vector[i][j] * distance_vector[i][j] *
                    inv_covariance_matrix[i][j][j];
    }
    /* do off diagonal once and count it twice */
    for_less(j, 0, num_features)
      for_less(k, j+1, num_features) {
        exponent += 2.0 * distance_vector[i][j] * 
                       inv_covariance_matrix[i][j][k] * distance_vector[i][k];
      } 
    
    fuzzy_bayes_vector[i] = normalize_class_vector[i] * exp(-0.5*exponent);

    if (apriori) { /* modify probability using spatial priors */

      /* fuzzy_bayes_vector[i] *= apriori_vector[i]; */
      fuzzy_bayes_vector[i] *= apriori_vector[mean_feature_class_vector[i]]; 
    }

    sum += fuzzy_bayes_vector[i];
  }

  /* normalize weights and select maximum */
  class_index = 0;
  max_prob = 0.0;
  if(sum > 0) {
    for_less(i, 0, num_classes) {
      fuzzy_bayes_vector[i] /= sum;
      /* check for new most probably class */
      if(fuzzy_bayes_vector[i] > max_prob) {
        max_prob = fuzzy_bayes_vector[i];
        class_index = i;
      }
    }
  }

  /* if fuzzy values are expected, copy them over */
  if ( class_prob )
    for_less( i, 0, num_classes ) {

      class_prob[i] = fuzzy_bayes_vector[i];
      class_labels[i] = mean_feature_class_vector[i]; 
    }


  /* set discrete classification */
  *class_num = mean_feature_class_vector[class_index];
}

/* ----------------------------- MNI Header -----------------------------------
@NAME       : bayesian_load_training
@INPUT      : 
@OUTPUT     : 
@RETURNS    : 
@DESCRIPTION: 
@METHOD     : 
@GLOBALS    : 
@CALLS      : 
@CREATED    : August 21, 1996 (John Sled)
@MODIFIED   : 
---------------------------------------------------------------------------- */
void bayesian_load_training(char *load_train_filename)
{
  int i, j, k;
  Real feature_value;

  FILE *learn_file;

  /* open the input training  file */
  learn_file = fopen(load_train_filename, "r");

  if ( learn_file == NULL) {
    fprintf(stderr, "Cannot open %s\n", load_train_filename);
    exit(EXIT_FAILURE);
  }

  /* scan for the number of features */
  fscanf( learn_file, "num_of_features = %d\n", &num_features);

  /* scan for the max number of classes M\n*/
  fscanf( learn_file, "num_of_classes = %d\n", &num_classes);

  if (debug) {
    fprintf( stderr, "num_of_features = %d\n", num_features);
    fprintf( stderr, "num_of_classes = %d\n", num_classes);
  }

  bayesian_allocate_memory();

  /* reserve a character pointer ( char *) for each class name */
  ALLOC( class_name, num_classes ); 
 
  /* load feature matrix from learn file */
  if (verbose)
    fprintf(stderr, "Loading the training file...\n");

  for_less( k, 0, num_classes) {
    
    fscanf(learn_file, "class = %d\n", &mean_feature_class_vector[k]);    

    /* allocate some space for each class_name[k], 5 for ex and
       simulate itoa, this could be replaced by a more elegant looking
       code, when I have time */
    ALLOC( class_name[k], 5);
    sprintf( class_name[k], "%d", mean_feature_class_vector[k]); 

    for_less( i, 0, num_features) {
      fscanf(learn_file, "%lf\n", &feature_value);
      mean_feature_matrix[k][i] = feature_value ;
    }

    for_less( i, 0, num_features) 
      for_less( j, 0, num_features) {
        fscanf(learn_file, "%lf\n", &feature_value);
        inv_covariance_matrix[k][i][j] = feature_value ;
    }

  }

  if (debug > 2 ) {

    fprintf( stderr, "Printing mean_feature_matrix ...\n");
    for_less( i, 0, num_classes) {
      for_less( j, 0, num_features) 
	fprintf( stderr, "%f ", mean_feature_matrix[i][j]);
      fprintf( stderr, "%d\n", mean_feature_class_vector[i]);      
    }

    fprintf( stderr, "-----\n");

    fprintf( stderr, "Printing inverse covariance_matrix ...\n");

    for_less( i, 0, num_classes) {
      fprintf( stderr, "\nCovariance for class %d\n", 
               mean_feature_class_vector[i]);
      for_less( j, 0, num_features) { 
        for_less( k, 0, num_features) 
          fprintf( stderr, "%lf ", inv_covariance_matrix[i][j][k]);
        fprintf(stderr, "\n");
      }
      fprintf( stderr, "-----\n");
    }
  }

  fclose(learn_file);

  /* put in scale factor for standard deviation */
  scale_matrices(inv_covariance_matrix, 
                 1.0/(stdev_scale_factor*stdev_scale_factor));
  if (debug > 2 ) {
    fprintf( stderr, "Scaling inverse covariance by %g\n",
             1.0/(stdev_scale_factor*stdev_scale_factor));
  }

  /* compute normalization factors */
  int l;
  for_less( l, 0, num_classes ) {
    normalize_class_vector[l] = 
      sqrt(matrix_determinant(num_features, inv_covariance_matrix[l]) /
           pow(2.0 * M_PI, num_features)); 
  }
}


/* ----------------------------- MNI Header -----------------------------------
@NAME       : bayesian_save_training
@INPUT      : 
@OUTPUT     : 
@RETURNS    : 
@DESCRIPTION: 
@METHOD     : 
@GLOBALS    : 
@CALLS      : 
@CREATED    : August 21, 1996 (John Sled)
@MODIFIED   : 
---------------------------------------------------------------------------- */
void bayesian_save_training(char *save_train_filename)
{

  int i, j, k;                 /* counters */
  FILE *train_file;         /* save filename */



  /* remove scale factor for standard deviation */
  scale_matrices(inv_covariance_matrix, 
                 (stdev_scale_factor*stdev_scale_factor));
  if (debug > 2 ) {
    fprintf( stderr, "Scaling inverse covariance by %g\n",
             1.0/(stdev_scale_factor*stdev_scale_factor));
  }

  train_file = fopen(save_train_filename, "w");


  /* see if the file could be opened */
  if ( train_file == NULL) {

    printf("Cannot open %s\n", save_train_filename);
    exit(EXIT_FAILURE);

  }

  /* store in the file the number of features trained on */
  fprintf(train_file, "num_of_features = %d\n", num_features);

  /* store in the file the max number of classes trained on */
  fprintf(train_file, "num_of_classes = %d\n", num_classes);
  
  for_less ( i, 0, num_classes ) {
    fprintf(train_file, "class = %d\n", mean_feature_class_vector[i]); 
    for_less ( j, 0, num_features )
      fprintf(train_file, "%lf ", mean_feature_matrix[i][j]); 
    fprintf(train_file, "\n\n");
    for_less ( j, 0, num_features )
      for_less( k, 0, num_features )
        fprintf(train_file, "%lf ", inv_covariance_matrix[i][j][k]); 
    fprintf(train_file, "\n");
  }

  fclose(train_file);


  /* put in scale factor for standard deviation */
  scale_matrices(inv_covariance_matrix, 
                 (stdev_scale_factor*stdev_scale_factor));
  if (debug > 2 ) {
    fprintf( stderr, "Scaling inverse covariance by %g\n",
             1.0/(stdev_scale_factor*stdev_scale_factor));
  }
}



/* ----------------------------- MNI Header -----------------------------------
@NAME       : matrix_determinant
@INPUT      : 
@OUTPUT     : 
@RETURNS    : 
@DESCRIPTION: 
@METHOD     : the simple slow method that you would do by hand
@GLOBALS    : 
@CALLS      : 
@CREATED    : August 21, 1996 (John Sled)
@MODIFIED   : 
---------------------------------------------------------------------------- */
Real matrix_determinant(int dimension, Real **matrix)
{
  Real determinant = 0;
  unsigned col;
  int factor, i, j;
  Real **sub_matrix;

  if (dimension > 1) {
    ALLOC2D(sub_matrix, dimension-1, dimension-1);
    for(col = 0, factor = 1; col < dimension; col++, factor *= -1) {
      /* create minor */
      for_less(i, 1, dimension) {
        for_less(j, 0, col) 
          sub_matrix[i-1][j] = matrix[i][j];
        for_less(j, col+1, dimension)
          sub_matrix[i-1][j-1] = matrix[i][j];
      }
      /* compute recursively */
      determinant += matrix[0][col] * 
        matrix_determinant(dimension-1, sub_matrix) * factor;
    }
    FREE2D(sub_matrix);
  }
  else 
    determinant = matrix[0][0];

  return determinant; 
}