lfcm.c

这是一个改进的快速实现模糊c-means聚类算法的程序

/* 
   lfcm.c
   A literal FCM implementation. 

   $Id: lfcm.c,v 1.3 2002/07/12 20:48:48 eschrich Exp $
   Steven Eschrich

   Copyright (C) 2002 University of South Florida

   This program is free software; you can redistribute it and/or modify it 
   under the terms of the GNU General Public License as published by the 
   Free Software Foundation; either version 2 of the License, or (at 
   your option) any later version.
   
   This program is distributed in the hope that it will be useful, but 
   WITHOUT ANY WARRANTY; without even the implied warranty of 
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 
   General Public License for more details.
   
   You should have received a copy of the GNU General Public License along 
   with this program; if not, write to the Free Software Foundation, Inc., 
   59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*/

#include <stdio.h>
#include <stdlib.h>
#include <sys/types.h>
#include <math.h>
#include <sys/times.h>
#include <sys/resource.h>
#include <limits.h>
#include <unistd.h>
#include <time.h>
#include <string.h>
#include "utils.h"





#define U(i,j) U[j][i]

/* Variables are defined statically here. These are reasonable
   defaults, but can be changed via the parameter to lfcmCluster() */
double epsilon=0.25;
double m=2.0;
int    C=2;
int    S=2;
int    N=2;
double **V;
double **U;
double **X;
int     number_of_iterations;
long seed;
int *max_value;



/* Public functions */
int lfcm();


/* Private functions */
int    update_centroids();
double update_umatrix();

/* Utilities */
int    init();
int    is_example_centroid();
double distance();

int output_centroids();
int output_umatrix();
int output_members();

/* External functions */
int load_test_data();
int load_atr_data();
int load_mri_data();


/* For testing purposes, we hard-code the desired number of clusters */
#define ATR_NUMBER_OF_CLUSTERS 5
#define MRI_NUMBER_OF_CLUSTERS 10
#define TEST_NUMBER_OF_CLUSTERS 2
#define TEST  1
#define ATR   2
#define MRI   3


/* Global variables */
int     dataset_type=MRI;
int     write_centroids=0;
int     write_umatrix=0;
int     write_members=0;

/* Variables that must be defined for called functions */
int  vals[][3]={{256,256,256},{0,0,0},{256,256,256},{4096,4096,4096}};


/* Function prototypes */
double *timing_of();  /* Calculate time in seconds */




int main(int argc, char **argv)
{

  struct rusage start_usage, end_usage;
  int ch;
  double *perf_times;
  char   *filename;

  epsilon=0.225;
  m=2.0;
  seed=2000;
  max_value=vals[dataset_type];


  while ( (ch=getopt(argc, argv,"hw:d:s:")) != EOF ) {
    switch (ch) {
    case 'h':
      fprintf(stderr,"Usage\n" \
	      "-d [a|t|m|s] Use dataset atr, mri, test, seawifs\n"\
	      "-w write cluster centers and memberships out\n"\
	      "-s seed  Use seed as the random seed\n");
      exit(1);
    case 'w':
      if ( !strcmp(optarg,"umatrix") ) write_umatrix=1;
      if ( !strcmp(optarg,"centroids") ) write_centroids=1;
      if ( !strcmp(optarg,"members") ) write_members=1;
      if ( !strcmp(optarg,"all")) 
	write_umatrix=write_centroids=write_members=1;
      break;
    case 'd':
      if ( *optarg == 'a' ) dataset_type=ATR;
      if ( *optarg == 'm' ) dataset_type=MRI;
      if ( *optarg == 't' ) dataset_type=TEST;
      max_value=vals[dataset_type];
      break;
    case 's':
      seed=atol(optarg);
      break;
    default:
    }
  }
  /* Print out main parameters for this run */
  fprintf(stdout,"FCM Parameters\n clusterMethod=literal fcm\n");
  filename=argv[optind];
  fprintf(stdout," file=%s\n\n",filename);



  /* Load the dataset, using one of a particular group of datasets. */
  switch (dataset_type) {
  case TEST:
    load_test_data(&X, &S, &N);
    C=TEST_NUMBER_OF_CLUSTERS;
    break;
  case ATR:
    load_atr_data(argv[optind],&X, &S, &N);
    C=ATR_NUMBER_OF_CLUSTERS;
    break;
  case MRI:
    load_mri_data(argv[optind], &X, &S, &N);
    C=MRI_NUMBER_OF_CLUSTERS;
    break;
  }


  fprintf(stdout, "Beginning to cluster...\n");


  /* Time the brfcm algorithm */
  getrusage(RUSAGE_SELF, &start_usage);
  lfcm();
  getrusage(RUSAGE_SELF, &end_usage);


  /* Output whatever clustering results we need */
  if ( write_centroids ) output_centroids(filename);
  if ( write_umatrix   ) output_umatrix(filename);
  if ( write_members   ) output_members(filename);

   
  /* Output timing numbers */
  perf_times=timing_of(start_usage, end_usage);
  printf("Timing: %f user, %f system, %f total.\n", 
				perf_times[0], perf_times[1], perf_times[0] +
				perf_times[1]);

  printf("Clustering required %d iterations.\n", number_of_iterations);




  return 0;
}






/* Main entry point into code. Cluster the dataset, given the details
   in the parameter block. */
int lfcm()
{
   double sqrerror = 2 * epsilon;

   /* Initialize code  */
   init();


   /* Run the updates iteratively */
   while (sqrerror > epsilon ) {
      number_of_iterations++;
      update_centroids();
      sqrerror=update_umatrix();
   }
   
   /* We go ahead and update the centroids - presumably this will not 
      change much, since the overall square error in U is small */
   update_centroids();

   return 0;
}

/* 
   update_centroids()
    Given a membership matrix U, recalculate the cluster centroids as the
    "weighted" mean of each contributing example from the dataset. Each
    example contributes by an amount proportional to the membership value.
*/
int update_centroids()
{
  int i,k,x;
  double numerator[S], denominator[S];

  /* For each cluster */
  for (i=0; i < C; i++)  {
   
    /* Zero out numerator and denominator options */
    for (x=0; x < S; x++) {
      numerator[x]=0;
      denominator[x]=0;
    }

    /* Calculate numerator */
    for (k=0; k < N; k++) {
      for (x=0; x < S; x++) 
	numerator[x] += pow(U(i,k), m) * X[k][x];
    }
    
    /* Calculate denominator */
    for (k=0; k < N; k++) {
      for (x=0; x < S; x++) 
	denominator[x] += pow(U(i,k), m);
    }
    
    /* Calculate V */
    for (x=0; x < S; x++) {
      V[i][x]= numerator[x] / denominator[x];
    }
    
  }  /* endfor: C clusters */

  return 0;
}

double update_umatrix()
{
  int i,j,k;
  int example_is_centroid;
  double summation, D_ki, D_kj;
  double square_difference=0;
  double newU;

  /* For each example in the dataset */
  for ( k=0; k < N; k++) {
    
    /* Special case: If Example is equal to a Cluster Centroid,
       then U=1.0 for that cluster and 0 for all others */
    if ( (example_is_centroid=is_example_centroid(k)) != -1 ) {
      for (i=0; i < C; i++) {
	if ( i == example_is_centroid )
	  U(i,k)=1.0;
	else 
	  U(i,k)=0.0;
      }
      continue;
    }
    
    /* For each class */
    for (i=0; i < C; i++) {
      summation=0;

      /* Calculate summation */
      for (j=0; j < C; j++) {
	D_ki=distance(X[k], V[i]);
	D_kj=distance(X[k], V[j]);
	summation += pow( D_ki / D_kj , (2.0/ (m-1)));
      }

      /* Weight is 1/sum */
      newU=1.0/(double)summation;
      
      /* Add to the squareDifference */
      square_difference += pow(U(i,k) - newU, 2);
      
      U(i,k)=newU;
    }

  } /* endfor n */
  

  return square_difference;

}

/*===================================================
  Utilities

  init()
  checkIfExampleIsCentroid()
  distance()

  ===================================================*/

/* Allocate storage for U and V dynamically. Also, copy over the
   variables that may have been externally set into short names,
   which are private and easier to access.
*/
int init()
{
  int i,j;
  
  /* Allocate necessary storage */
  V=(double **)CALLOC(C,sizeof(double *));
  for (i=0; i < C; i++) 
    V[i]=(double *)CALLOC(S, sizeof(double));
  
  U=(double **)CALLOC(N, sizeof(double *));
  for (i=0; i < N; i++)
    U[i]=(double *)CALLOC(C,sizeof(double));


  /* Place random values in V, then update U matrix based on it */
  srand48(seed);
  for (i=0; i < C; i++) {
    for (j=0; j < S; j++) {
      V[i][j]=drand48() * max_value[j];
    }
  }
  /* Once values are populated in V, update the U Matrix for sane values */
  update_umatrix();

  return 0;
}

  
/* If X[k] == V[i] for some i, then return that i. Otherwise, return -1 */
int is_example_centroid(int k)
{
  int  i,x;

  for (i=0; i < C; i++) {
    for (x=0; x < S; x++) {
      if ( X[k][x] != V[i][x] ) break;
    }
    if ( x == S )  /* X==V */
      return i;
  }
  return -1;
}

double distance(double *v1, double *v2)
{
  int x;
  double sum=0;

  for (x=0; x < S; x++) 
    sum += (v1[x] - v2[x]) * (v1[x] - v2[x]);

  return sqrt(sum);
}



/*=====================================================
  Public output utilities 

  output_centroids()
  output_umatrix()
  output_members()
  =====================================================*/
int output_centroids(char *filestem)
{
  FILE *fp;
  char buf[1024];
  int i,j;
  
  sprintf(buf,"%s.centroids", filestem);
  fp=FOPEN(buf,"w");
  for (i=0;i < C ;i++) {
    for (j=0; j < S; j++)
       fprintf(fp, "%f\t",V[i][j]);
    fprintf(fp,"\n");
  }
  fclose(fp);

  return 0;
}

int output_umatrix(char *filestem)
{
  FILE *fp;
  char buf[1024];
  int i,j;

  sprintf(buf,"%s.umatrix", filestem);
  fp=FOPEN(buf,"w");
  for (i=0; i < N; i++) {
    for (j=0; j < C; j++)
      fprintf(fp,"%f\t", U[i][j]);
    fprintf(fp,"\n");
  }
  fclose(fp);

  return 0;
}

int output_members(char *filestem)
{
  FILE *fp;
  char buf[1024];
  int i,j,max;

  sprintf(buf,"%s.members", filestem);
  fp=FOPEN(buf,"w");
  for (i=0; i < N; i++) {
    for (max=j=0; j < C; j++)
      if ( U[i][j] > U[i][max] ) max=j;
    fprintf(fp,"%d\n",max);
  }
  fclose(fp);

  return 0;
}