cluster1.c

it is the Data Mining Algorithm source code.

/*----------------------------------------------------------------------
  File    : cluster1.c
  Contents: cluster and cluster set management (basic functions)
  Author  : Christian Borgelt
  History : 2001.09.05 file created
            2001.09.11 functions cls_range, cls_reg, cls_expand added
            2001.09.15 first version of fuzzy c-means completed
            2002.08.07 parsing function completed
            2003.01.30 data normalization code added
            2003.02.01 extended functions (for tuples) added
            2003.02.02 extended parse function added
            2003.02.07 reinitialization added to function cls_range
            2003.03.09 some initialization moved to separate routines
            2003.03.13 membership degree normalization made optional
            2003.04.10 matrix initialization bug in _params fixed
            2003.05.15 membership normalization modes added
            2003.05.16 noise clustering added
            2003.08.06 bug in mapping binary attributes fixed
            2003.08.11 adapted to new module attmap
            2003.08.12 adapted to new module nstats
            2003.08.14 bug in _params fixed (covariance matrix reading)
            2003.08.15 adapted to new module radfn
            2004.02.16 bug in _params for array of variances fixed
            2004.02.23 bug in function _params fixed
            2004.02.25 adapted to changed radfn, bug in cls_exec fixed
            2004.02.27 function cls_deletex added
            2004.03.01 update functions moved to a separate file
            2004.03.17 missing evaluation of cluster weight added
            2004.03.18 normalization parameters added
            2004.04.25 bug in function cls_unscale fixed
            2004.06.21 reading covariances made more robust
            2004.07.28 bug in function _params (cls_parse) fixed
            2004.08.12 adapted to new module parse
            2007.01.10 use of attribute maps changed, bug in _norm fixed
            2007.11.26 noise cluster configured with separate function
            2007.11.27 membership degree transformation extended
            2008.02.29 redesign for feature weighting completed
            2008.04.02 reading and writing of joint (co)variances added
            2008.04.03 numerical problems solved, consistent matrices
            2008.04.08 auxiliary functions extended (cls_vars etc.)
            2008.04.11 rescaling removed from covariance matrices
            2008.04.14 cls_unscale adapted to removed rescaling
            2008.04.25 check for positive eigenvalues added to _covmat
----------------------------------------------------------------------*/
#include <stdlib.h>
#include <limits.h>
#include <float.h>
#include <string.h>
#include <assert.h>
#include "cluster.h"
#include "vecops.h"

/*----------------------------------------------------------------------
  Preprocessor Definitions
----------------------------------------------------------------------*/
#define BLKSIZE     16          /* block size for cluster array */
#define MINVAR      1e-48       /* minimal variance */
#define MAXVAR      1e+48       /* maximal variance */
#define MINDET      1e-48       /* minimal determinant */
#define MAXDET      1e+48       /* maximal determinant */

/*----------------------------------------------------------------------
  Constants
----------------------------------------------------------------------*/
static char *nrmnames[] =       /* names of normalization modes */
  { "none", "sum1", "max1", "hard" };

/*----------------------------------------------------------------------
  Auxiliary Functions
----------------------------------------------------------------------*/

static CLSET* _create (int incnt, int clscnt)
{                               /* --- create the base structure */
  CLSET *clset = (CLSET*)malloc(sizeof(CLSET));
  if (!clset) return NULL;      /* create the base structure */
  clset->type     = CLS_CENTER; /* prototypes consist of centers only */
  clset->tnew     = CLS_CENTER;
  clset->incnt    = incnt;
  clset->clscnt   = clscnt;
  clset->init     = 0;
  clset->cls      = NULL;
  clset->radfn    = rf_cauchy;  /* default radial function: */
  clset->drvfn    = rfd_cauchy; /* \mu(x) = 1 / d^2(\vec{x}, \vec{c}) */
  clset->rfnps[0] = 2;          /* (standard fuzzy clustering) */
  clset->rfnps[1] = 0;
  clset->unit     = 1;
  clset->fwexp    = 0;          /* no feature weighting/selection */
  clset->noise    = 0;          /* do not use a noise cluster */
  clset->ncmsd[0] = 0;
  clset->ncmsd[1] = 0;
  clset->norm     = CLS_SUM1;   /* default normalization: */
  clset->nrmps[0] = 1;          /* u_ij = u_ij^* / sum_k u_{kj}^* */
  clset->nrmps[1] = 0;
  clset->msexp    = 2;          /* use standard fuzzy clustering */
  clset->method   = CLS_ALTOPT; /* with alternating optimization */
  clset->moment   = 0;
  clset->growth   = 1.2;
  clset->shrink   = 0.5;
  clset->maxchg   = 1;
  clset->maxfac   = clset->growth /(clset->growth+1);
  clset->rates[0] = 0.2;
  clset->rates[1] = 0.2;
  clset->rates[2] = 0.2;
  clset->decay[0] = -0.2;
  clset->decay[1] = -0.2;
  clset->decay[2] = -0.2;
  clset->regps[0] = 0;          /* regularization parameters */
  clset->regps[1] = 1;          /* (b, a, s, h, w) */
  clset->regps[2] = 1;
  clset->regps[3] = 0;
  clset->regps[4] = 0;
  clset->steps    = 0;
  clset->resized  = 0;
  clset->nst      = NULL;       /* initialize pointers, */
  clset->vec      = NULL;       /* so that deletion with */
  clset->buf      = NULL;       /* cls_delete works properly */
  clset->mat      = NULL;
  return clset;                 /* return the created cluster set */
}  /* _create() */

/*--------------------------------------------------------------------*/

static int _initcls (CLUSTER *c, int dim)
{                               /* --- init. a cluster prototype */
  assert(c && (dim > 0));       /* check the function arguments */
  c->cov = mat_create(dim, dim);
  c->inv = mat_create(dim, dim);
  c->smp = mat_create(dim, dim);
  c->chv = mat_create(dim, dim);
  c->grv = mat_create(dim, dim);
  if (!c->cov || !c->inv || !c->smp || !c->chv || !c->grv)
    return -1;                  /* create cluster-specific matrices */
  c->ctr = mat_vector(c->cov);  /* get their associated vectors */
  c->dif = mat_vector(c->inv);  /* for easier access */
  c->sum = mat_vector(c->smp);
  c->chc = mat_vector(c->chv);
  c->grc = mat_vector(c->grv);
  return c->dec = 0;            /* return 'ok' */
}  /* _initcls() */

/*--------------------------------------------------------------------*/

void cls_combine (CLSET *clset)
{                               /* --- combine (co)variances */
  int     i, m;                 /* loop variable, matrix location */
  CLUSTER *c, *r;               /* to traverse the clusters */
  double  w;                    /* cluster weight */

  assert(clset                  /* check the function argument */
     && (clset->type & CLS_JOINT));
  if      (clset->type & CLS_COVARS) m = MAT_UPPER;
  else if (clset->type & CLS_VARS)   m = MAT_DIAG;
  else if (clset->type & CLS_SIZE)   m = MAT_CORNER;
  else return;                  /* get the matrix location */
  r = clset->cls;               /* get the reference cluster */
  mat_mulsx(r->cov, r->cov, mat_getwgt(r->cov), m);
  for (c = r +(i = clset->clscnt); --i > 0; ) {
    w = mat_getwgt((--c)->cov); /* traverse the clusters */
    mat_addx(r->cov, r->cov, w, c->cov, m);
  }                             /* sum the (co)variances weighted */
  /* It is assumed that the sum of the cluster weights is 1 and */
  /* thus that no renormalization of the result is necessary.   */
}  /* cls_combine() */

/*--------------------------------------------------------------------*/

void cls_distrib (CLSET *clset)
{                               /* --- distribute (co)variances */
  int     i, m;                 /* loop variable, matrix location */
  CLUSTER *c, *r;               /* to traverse the clusters */

  assert(clset                  /* check the function argument */
     && (clset->type & CLS_JOINT));
  if      (clset->type & CLS_COVARS)          m = MAT_FULL;
  else if (clset->type & (CLS_VARS|CLS_SIZE)) m = MAT_DIAG;
  else return;                  /* get the matrix location */
  r = clset->cls;               /* get the reference cluster */
  for (c = r +(i = clset->clscnt); --i > 0; ) {
    (--c)->var = r->var;        /* traverse the clusters and */
    c->size    = r->size;       /* copy the isotropic variance, */
    c->scale   = r->scale;      /* the cluster size and the scaling */
    mat_copyx(c->cov, r->cov, m);
    mat_copyx(c->inv, r->inv, m);
  }                             /* copy (co)variances and inverses */
}  /* cls_distrib() */

/*--------------------------------------------------------------------*/

static void _covmat (CLUSTER *c, int eigen)
{                               /* --- decompose covariance matrix */
  int    i, k, n;               /* loop variable, number of dims. */
  double t, det;                /* buffers for shift, determinant */

  assert(c);                    /* check the function arguments */
  n = mat_rowcnt(c->cov);       /* get the number of dimensions */
  for (t = 0.5*MINVAR, i = 96; --i >= 0; ) {
    if (c->dec) k = 0;          /* do several decomposition trials */
    else {                      /* if decomposition does not exist */
      if (!eigen)               /* do Cholesky decomposition */
        k = mat_chdecom(c->inv, c->cov);
      else {                    /* do eigenvalue decomposition */
        mat_3dred(c->buf, c->buf +n, c->inv, c->cov, MAT_UPPER);
        k = mat_3dqli(c->dif, c->inv, c->buf, c->buf +n, n, 256);
        if (k == 0) {           /* if successfully decomposed */
          mat_evsort(c->dif, c->inv, 1);
          k = (c->dif[0] >= MINVAR) ? 0 : -1;
        }                       /* check for positive eigenvalues */
      }                         /* (positive definite matrix) */
    }
    if (k == 0) {               /* if decomposition was successful */
      if (eigen) {              /* if eigenvector decomposition */
        det    = vec_prod(c->dif, n);
        c->var = ((det >= MINDET) && (det <= MAXDET))
               ? pow(det, 1.0/n) : exp(vec_sumlog(c->dif, n) /n); }
      else {                    /* if Cholesky decomposition */
        det    = mat_chdet(c->inv);
        c->var = ((det >= MINDET) && (det <= MAXDET))
               ? pow(det, 1.0/n) : exp(mat_chlogd(c->inv) /n);
      }                         /* compute the isotropic variance */
      if (c->var > MINVAR) break;
    }                           /* check for a proper value */
    mat_diaadds(c->cov, t += t);/* if the decomposition failed, */
  }                             /* shift the eigenvalues */
  if (i < 0) {                  /* if decomposition failed totally, */
    c->var = 1;                 /* set a unit covariance matrix */
    mat_init(c->cov, MAT_UNIT, NULL);
    mat_init(c->inv, MAT_UNIT, NULL);
    if (eigen) { while (--n >= 0) c->dif[i] = 1; }
  }                             /* set unit eigenvalues if necessary */
  c->dec = -1;                  /* set the decomposition flag */
}  /* _covmat() */

/*--------------------------------------------------------------------*/

static void _diamat (CLUSTER *c)
{                               /* --- diagonal covariance matrix */
  int    i, n;                  /* loop variable, number of dims. */
  double t, det;                /* buffers for shift, determinant */

  assert(c);                    /* check the function argument */
  n = mat_rowcnt(c->cov);       /* get the number of dimensions */
  for (t = MINVAR, i = 97; --i >= 0; ) {
    det = mat_diaprod(c->cov);  /* compute the determinant */
    c->var = ((det >= MINDET) && (det <= MAXDET))
           ? pow(det, 1.0/n) : exp(mat_dialog(c->cov) /n);
    if (c->var > MINVAR) break; /* compute the isotropic variance */
    mat_diaadds(c->cov, t);     /* if the computation failed, */
    if (i < 96) t += t;         /* shift the eigenvalues and */
  }                             /* double the increment value */
  if (i >= 0) return;           /* if sucessful computation, abort */
  c->var = 1;                   /* otherwise set a unit cluster size */
  mat_diasetx(c->cov, 1);       /* and a unit covariance matrix */
}  /* _diamat() */

/*----------------------------------------------------------------------