cluster-sim.c

经典cure聚类算法,实现用c++语言.大家

/*
 * cluster.c -- cluster main program
 *
 * $Log: cluster.c,v $
 * Revision 1.23  1993/02/03  07:43:07  stolcke
 * added code vector output for cluster
 *
 * Revision 1.22  1993/01/20  10:28:02  stolcke
 * faster O(n^2) clustering
 *
 * Revision 1.21  1993/01/16  01:44:34  stolcke
 * use triangular distance matrix (save 1/2 the memory)
 *
 * Revision 1.20  1993/01/16  00:27:07  stolcke
 * make dflag printing more efficient
 *
 * Revision 1.19  1993/01/16  00:03:30  stolcke
 * avoid unnecessary reallocation of patterns
 * skip non-root distance entries more efficiently
 *
 * Revision 1.18  1992/02/12  07:13:13  stolcke
 * -E option added
 *
 * Revision 1.17  91/07/15  12:56:00  stolcke
 * width option added 
 * 
 * Revision 1.16  91/07/14  01:10:54  stolcke
 * curses support added
 * graphing routines moved to separate file
 * 
 * Revision 1.15  91/07/10  21:26:54  stolcke
 * saber-cleaned and bad bug in allocation macro fixed
 * 
 * Revision 1.14  91/04/20  16:17:59  stolcke
 * second release (beta)
 * 
 * Revision 1.13  91/04/18  19:03:39  stolcke
 * some more error checks
 * 
 * Revision 1.12  91/04/18  18:28:21  stolcke
 * general cleanup and partial rewrite
 * 
 * Revision 1.11  91/04/18  13:29:29  stolcke
 * merged pca into cluster
 *
 * Revision 1.10  91/04/12  02:09:48  stolcke
 * name quoting added
 *
 * Revision 1.9  90/12/06  21:47:46  stolcke
 * fixed averaging with D/C values
 *
 * Revision 1.8  90/12/06  14:32:52  stolcke
 * support for don't care values added
 * fixed bug when total_distance = 0
 *
 * Revision 1.7  90/11/12  20:31:22  stolcke
 * two initialization fixed
 *
 * Revision 1.6  90/11/12  18:43:26  stolcke
 * workaround to avoid old Sun compiler bug
 *
 * Revision 1.5  90/11/12  16:40:34  stolcke
 * added metric selection and scaling
 *
 * Revision 1.4  90/11/11  18:14:27  stolcke
 * -n flag no longer needed
 *
 * Revision 1.3  90/11/11  02:30:14  stolcke
 * no more seek needed, can use in pipe now
 *
 * Revision 1.2  90/11/10  23:44:02  stolcke
 * added new features (see manpage)
 *
 */

#if !defined(lint) && !defined(SABER)
static char rcsid[] = "$Header: /usr/src/local/conn/cluster/RCS/cluster.c,v 1.23 1993/02/03 07:43:07 stolcke Exp $";
#endif				/* not lint */

#include <stdio.h>
#include <string.h>
#include <math.h>
#include "alloc.h"
#include "error.h"
#include "cluster.h"

#define NONE (-2)

#define BUFSIZE 256

#ifndef SCALE
#define SCALE "_SCALE_"
#endif

#ifndef DONTCARE
#define DONTCARE "D/C"
#endif

static FLOAT   distance();
static FLOAT   similarity();
static FLOAT   CombSim();
static FLOAT   root();

static char buffer[BUFSIZE];	/* temporary string read buffer */

#define PCA	"pca"

/* different sim combining method */
#define SINGLE 0
#define COMPLETE 1
#define UPGMA 2
#define WPGMA 3
#define UPGMC 4
#define WPGMC 5
#define WARD 6

char   *program;		/* program name */

int     vflag = 0;		/* print explanatory messages */
static int pflag = 0;		/* do PCA instead of HCA */
static int dflag = 0;		/* print distances */
static int tflag = 0;		/* print tree in ASCII */
static int Tflag = 0;		/* display tree using curses */
static int gflag = 0;		/* print tree in graph(1) format */
static int bflag = 0;		/* print tree in graph(1) format, with breaks */
static int Bflag = 0;		/* print patterns as bit vectors */
static int sflag = 0;		/* suppress scaling  */
int	Eflag = 0;		/* print eigenvalues */

static int width = 0;		/* width of ASCII tree representation */

static int norm = 2;		/* which l-norm to use for distance metric */

static int csize = 0; 		/* number of clusters */

static 
usage()
{
    if (!pflag)
	fprintf(stderr,
	   "usage: %s [-dtTgbBsv] [-w width] [-n norm] [-k no-cluster] [vectorfile [namesfile]]\n", program);
    else
	fprintf(stderr,
		"usage: %s [-Esv] [-e eigenbase] [-c pcs] [vectorfile [namesfile]]\n", PCA);
    exit(2);
}

main(argc, argv)
    int     argc;
    char   *argv[];
{
    FILE   *fp, *pfile;
    char   *efile = NULL;
    char   *comps = NULL;
    FLOAT **pattern = NULL;
    int     lpat, npat;
    char  **name = NULL;
    int   i, *pinfo;
    char  fname[80];

    int     opt;
    extern char *optarg;
    extern int optind;

    /* who are we ? */
    if ((program = strrchr(argv[0], '/')) == NULL)
	program = argv[0];
    else
	program += 1;

    if (strcmp(program, PCA) == 0)
	pflag = 1;

    while ((opt = getopt(argc, argv, "pdtTgbBvsnk:e:c:w:E")) != -1)
	switch (opt) {
	case 'p':
	    pflag = 1;
	    break;
	case 'd':
	    if (pflag)
		usage();
	    dflag = 1;
	    break;
	case 't':
	    if (pflag)
		usage();
	    tflag = 1;
	    break;
	case 'T':
	    if (pflag)
		usage();
	    Tflag = 1;
	    break;
	case 'g':
	    if (pflag)
		usage();
	    gflag = 1;
	    break;
	case 'b':
	    if (pflag)
		usage();
	    bflag = 1;
	    gflag = 1;
	    break;
	case 'B':
	    if (pflag)
		usage();
	    Bflag = 1;
	    break;
	case 'v':
	    vflag = 1;
	    break;
	case 'n':
	    if (pflag)
		usage();
	    if (sscanf(optarg, "%d", &norm) != 1 || norm < 0)
		usage();
	    break;
	case 'k':
	    if (pflag)
		usage();
	    if (sscanf(optarg, "%d", &csize) != 1 || csize < 0)
		usage();
	    break;
	case 'e':
	    if (!pflag)
		usage();
	    efile = optarg;
	    break;
	case 'c':
	    if (!pflag)
		usage();
	    comps = optarg;
	    break;
	case 's':
	    sflag = 1;
	    break;
	case 'w':
	    if (pflag)
		usage();
	    if (sscanf(optarg, "%d", &width) != 1 || width <= 0)
		usage();
	    break;
	case 'E':
	    if (!pflag)
		usage();
	    Eflag = 1;
	    break;
	case '?':
	    usage();
	    break;
	}

    if (!(pflag || dflag || tflag || Tflag || gflag || Bflag))
	dflag = tflag = vflag = 1;	/* default behavior */

    if (optind + 2 < argc)
	usage();

    if (!(optind < argc) || !strcmp(argv[optind], "-"))
	fp = stdin;
    else
	IfErr(fp = fopen(argv[optind], "r")) {
	    fprintf(stderr, "%s: cannot open file %s\n", argv[0], argv[optind]);
	    exit(1);
	}

    IfErr(read_pattern(fp, &pattern, &lpat, &npat, &name)) {
	fprintf(stderr, "%s: %s: cannot read pattern\n", program, ERR_MSG);
	exit(1);
    }
    if (vflag)
	fprintf(stderr, "read %d patterns:  size = %d\n", npat, lpat);

    if (optind + 1 < argc) {
	IfErr (name = new_array_of(npat, char *)) {;
	    fprintf(stderr, "%s: not enough core for name array\n", program);
	    exit(1);
	}

	if (!strcmp(argv[optind + 1], "-"))
	    fp = stdin;
	else
	    IfErr(fp = fopen(argv[optind + 1], "r")) {
		fprintf(stderr, "%s: cannot open file %s\n",
			program, argv[optind + 1]);
		exit(1);
	    }

	IfErr(read_names(fp, name, npat)) {
	    fprintf(stderr, "%s: %s: cannot read names\n", program, ERR_MSG);
	    exit(1);
	}
    }

    IfErr(pinfo = new_array_of(npat, int)) {
	    fprintf(stderr, "%s: not enough core\n", program);
	    exit(1);
    }

    IfErr(pinfo = new_array_of(npat, int)) {
	    fprintf(stderr, "%s: not enough core\n", program);
	    exit(1);
    }

    sprintf(fname, "%s-partition", argv[optind]);
    IfErr(pfile = fopen(fname, "w")) {
	    fprintf(stderr, "%s: cannot open file %s\n", argv[0], argv[optind]);
	    exit(1);
    }

    cluster(pattern, name, lpat, npat, pinfo);

    for (i=0; i<npat; ++i)
	fprintf(pfile, "%d\n", pinfo[i]);
    free(pinfo);
    fclose(pfile);

    exit(0);
}

/* skip blanks and next end of line */
skip_blanks(fp)
    FILE   *fp;
{
    char    c;

    while ((c = getc(fp)) == ' ' || c == '\t');
    if (c != '\n')
	ungetc(c, fp);
    return c;
}

read_names(fp, name, npat)
    FILE   *fp;
    char  **name;
    int     npat;
{
    register int i;

    for (i = 0; i < npat; i++) {
	char    c = skip_blanks(fp);

	if (c == '\"') {
	    getc(fp);
	    fgets(buffer, sizeof(buffer), fp);
	    buffer[strlen(buffer) - 1] = '\0';
	}
	else {
	    IfEOF(fscanf(fp, "%s", buffer))
		Erreturn("not enough names");
	    skip_blanks(fp);
	}

	IfErr(name[i] = new_string(buffer))
	    Erreturn("not enough core");
    }
    return MY_OK;
}

void
print_names(tree, name)
    BiTree *tree;
    char  **name;
{
    if (tree->leaf != LEAF) {
	if (name)
	    printf(" %s", name[tree->leaf]);
	else
	    printf(" %d", tree->leaf);
    }
    else {
	print_names(tree->r_tree, name);
	print_names(tree->l_tree, name);
    }
}

void
print_k_cluster(tree, name, pinfo, cno)
    BiTree *tree;
    char  **name;
    int   *pinfo;
    int   cno;
{
    if (tree->leaf != LEAF) {
	if (name)
	    printf("%d : %s\n", tree->leaf, name[tree->leaf]);
	else
	    printf("%d : %d\n", tree->leaf, tree->leaf);

	pinfo[tree->leaf] = cno;
    }
    else {
	print_k_cluster(tree->r_tree, name, pinfo, cno);
	print_k_cluster(tree->l_tree, name, pinfo, cno);
    }
}


find_nnb(items, which, lpat, index, nsim)
    BiTree *items;	/* node array */
    int which;		/* index of node to find nearest neightbor for */
    int lpat;		/* pattern length */
    int *index;		/* returns: index of nnb (-1 if none) */
    FLOAT *nsim;	/* returns: sim to nnb */
{
    int i;
    FLOAT sim, max_sim;
    int max_index;

    if (items[which].root == FALSE) {
	*index = NONE;
	return;
    }

    max_index = NONE;
    max_sim = 0.0;
    /*
     * find minimum distance neighbor -- to avoid duplication
     * only pairs with 1st index < 2nd index are considered
     */
    for (i = 0; i < which; i++) {
	if (items[i].root == FALSE)
	    continue;
	
	sim = items[which].slist[i];
	if (max_index == NONE || sim > max_sim) {
	    max_index = i;
	    max_sim = sim;
	}
    }

    if (max_index >= 0) {
        *index = max_index;
	*nsim = max_sim;
    }
    else {
	*index = NONE;
	*nsim = NONE;
    }
    return;
}

cluster(pattern, name, lpat, npat, pinfo)
    FLOAT **pattern;
    char  **name;
    int     lpat, npat;
    int   *pinfo;
{
    register int i, j, k, nodes = npat, cno;
    BiTree *item = new_array_of(npat, BiTree);
    /*
     * for each data point or cluster center, we keep the index of the nearest
     * neighbor, as well as the distance to it.
     */
    int *nnb_index = new_array_of(npat, int);