cluster2.c

it is the Data Mining Algorithm source code.

/*----------------------------------------------------------------------
  File    : cluster2.c
  Contents: cluster and cluster set management (update functions)
  Author  : Christian Borgelt
  History : 2001.09.05 file created as cluster1.c
            2001.09.12 function cls_init completed
            2001.09.15 first version of fuzzy c-means completed
            2001.09.16 hard c-means clustering algorithm added
            2002.09.09 neural network update methods added
            2003.01.29 some cleanup of neural network methods
            2003.01.31 initialization for resilient method changed
            2003.02.18 multiple init. in mode CLS_POINTS removed
            2003.06.07 cluster size adaptation added
            2003.10.29 bug in function cls_init fixed (CLS_POINTS)
            2004.02.22 size computation for spherical clusters modified
            2004.02.23 bug in function _quick fixed
            2004.02.29 some bugs in cluster size computation fixed
            2004.03.01 update functions moved to this file
            2004.03.02 shape and size regularization added
            2004.03.19 weight/prior regularization added
            2004.04.12 competitive learning function completed
            2004.04.13 regularization adapted for competitive learning
            2004.04.14 treatment of (almost) empty clusters improved
            2004.04.15 update function improved (loop moved)
            2004.04.23 some numeric problems with high-dim. data solved
            2004.04.26 more numeric problems with high-dim. data solved
            2004.04.27 upper learning rate bound added (alt. est. step)
            2004.04.30 rescaling removed from parameter est. functions
            2004.07.13 normalization of center vectors added
            2004.07.15 bug in center normalization for _complrn fixed
            2004.07.28 update of centers only added to cls_update
            2004.08.14 bug in initialization with CLS_POINTS fixed
            2004.08.18 first version of backpropagation functions added
            2007.11.22 eigenvalue sorting added for covariance mode
            2007.11.27 unit centers and centers at origin corrected
            2008.02.29 redesign for feature weighting completed
            2008.04.08 regularization functions adapted to new version
            2008.04.11 neural network inspired update debugged
            2008.04.14 neural network inspired update improved
            2008.04.21 modify only covariance changes towards zero
            2008.04.22 momentum term update bounded by maximal change
            2008.04.25 covariance matrix checks improved in _neural
----------------------------------------------------------------------*/
#include <stdlib.h>
#include <float.h>
#include <math.h>
#include <assert.h>
#include "cluster.h"
#include "vecops.h"

/*----------------------------------------------------------------------
  Preprocessor Definitions
----------------------------------------------------------------------*/
#define MINVAR      1e-48       /* minimal variance */
#define MINDET      1e-48       /* minimal determinant */
#define MAXDET      1e+48       /* maximal determinant */
#define MINWEIGHT   1e-12       /* minimal cluster weight */

/*----------------------------------------------------------------------
  Type Definitions
----------------------------------------------------------------------*/
typedef MATRIX* MATADDFN (MATRIX *mat, const double *vec, double wgt);
typedef double  UPDATEFN (CLSET *clset, double  grd,
                          double prv,   double *chg);

/*----------------------------------------------------------------------
  Auxiliary Functions from cluster1.c
----------------------------------------------------------------------*/
extern void cls_combine (CLSET *clset);
extern void cls_distrib (CLSET *clset);
extern void cls_vars    (CLSET *clset, int eigen);
extern void cls_gauge   (CLSET *clset);
extern void cls_resize  (CLSET *clset, int decomp);
extern void cls_ftwgts  (CLSET *clset);
extern void cls_invert  (CLSET *clset);
extern void cls_reinit  (CLSET *clset);

/*----------------------------------------------------------------------
  Gradient Function
----------------------------------------------------------------------*/

static void _gradient (CLSET *clset)
{                               /* --- gradient based update */
  /* ... to be done ... */
}  /* _gradient() */

/*----------------------------------------------------------------------
  Alternating Optimization Function
----------------------------------------------------------------------*/

static void _altopt (CLSET *clset)
{                               /* --- alternating optimization */
  int     i;                    /* loop variable */
  CLUSTER *c;                   /* to traverse the clusters */
  MATRIX  *mat;                 /* exchange buffer for matrices */
  double  w, sum = 0;           /* (sum of) cluster weights */

  assert(clset);                /* check the function argument */

  /* --- check and sum the weights --- */
  for (c = clset->cls +(i = clset->clscnt); --i >= 0; ) {
    w = mat_getwgt((--c)->smp); /* traverse the clusters */
    if (w >= MINWEIGHT) {                     c->ok = -1; }
    else { mat_setwgt(c->smp, w = MINWEIGHT); c->ok =  0; }
    sum += w;                   /* check and sum the cluster weights */
    mat = c->smp; c->smp = c->cov; c->cov = mat;
    c->ctr = mat_vector(c->cov);/* exchange the covariance matrix */
    c->sum = mat_vector(c->smp);/* with the aggregation matrix and */
  }                             /* get the corresponding buffers */
  sum = 1/sum;                  /* compute the normalization factor */

  /* --- compute new centers and (co)variances --- */
  c += i = clset->clscnt;       /* traverse the clusters */
  if      (clset->type & CLS_COVARS) {
    while (--i >= 0) { --c;     /* if adaptable covariances */
      if (c->ok) mat_covarx(c->cov, c->cov, 1);
      else mat_copyx(c->cov, c->smp, MAT_VECTOR|MAT_UPPER);
      mat_mulwgt(c->cov, sum);  /* compute new covariance matrix */
    } }                         /* and normalize the weights */
  else if (clset->type & CLS_VARS) {
    while (--i >= 0) { --c;     /* if adaptable variances */
      if (c->ok) mat_varx(c->cov, c->cov, 1);
      else mat_copyx(c->cov, c->smp, MAT_VECTOR|MAT_DIAG);
      mat_mulwgt(c->cov, sum);  /* compute new variances */
    } }                         /* and normalize the weights */
  else if (clset->type & CLS_SIZE) {
    while (--i >= 0) { --c;     /* if adaptable size */
      if (c->ok) mat_isovarx(c->cov, c->cov, 1);
      else mat_copyx(c->cov, c->smp, MAT_VECTOR|MAT_CORNER);
      mat_mulwgt(c->cov, sum);  /* compute new cluster size */
    } }                         /* and normalize the weight */
  else {                        /* if no shape or size flags set */
    while (--i >= 0) { --c;     /* (fixed isotropic variance) */
      if (c->ok) mat_mean(c->ctr, c->cov);
      else mat_copyx(c->cov, c->smp, MAT_VECTOR);
      mat_mulwgt(c->cov, sum);  /* only compute new cluster centers */
    }                           /* and normalize the weights */
  }                             /* (copy old values if necessary) */

  /* --- combine (co)variances --- */
  if ((clset->type &  CLS_JOINT)
  &&  (clset->type & (CLS_COVARS|CLS_VARS|CLS_SIZE)))
    cls_combine(clset);         /* combine the (co)variances */
}  /* _altopt() */

/*----------------------------------------------------------------------
  Competitive Learning Function
----------------------------------------------------------------------*/

static void _complrn (CLSET *clset)
{                               /* --- competitive learning */
  int     i, m;                 /* loop variable, matrix mode */
  CLUSTER *c;                   /* to traverse the clusters */
  double  w, sum;               /* (sum of) cluster weights */
  double  lrn, dec, eta;        /* learning rate and related values */
  MATRIX  *mat;                 /* exchange buffer */

  assert(clset);                /* check the function argument */

  /* --- compute new centers --- */
  dec = clset->decay[0];        /* get the decay parameter and */
  lrn = clset->rates[0];        /* compute the next learning rate */
  if ((dec > 0) && (dec < 1)) lrn *= pow(dec, clset->steps);
  else if          (dec < 0)  lrn *= pow(clset->steps+1, dec);
  for (c = clset->cls +(i = clset->clscnt); --i >= 0; ) {
    w = mat_getwgt((--c)->smp); /* traverse the clusters */
    w = (w > 1) ? lrn/w : lrn;  /* compute the learning rate */
    mat = c->smp; c->smp = c->cov; c->cov = mat;
    c->ctr = mat_vector(c->cov);/* exchange the covariance matrix */
    c->sum = mat_vector(c->smp);/* with the aggregation matrix */
    vec_add(c->ctr, clset->incnt, c->sum, w, c->ctr);
  }                             /* compute new center coordinates */
  /* Note that c->sum is the aggregate of the difference vectors to  */
  /* the cluster centers and *not* the aggregate of the data vectors */
  /* as for alternating estimation/fuzzy clustering. Hence there is  */
  /* no decay factor for the old cluster center c->ctr (cf. below).  */

  /* --- compute new (co)variances --- */
  if (clset->type & (CLS_COVARS|CLS_VARS|CLS_SIZE)) {
    dec = clset->decay[1];      /* get the decay parameter and */
    lrn = clset->rates[1];      /* compute the next learning rate */
    if ((dec > 0) && (dec < 1)) lrn *= pow(dec, clset->steps);
    else if          (dec < 0)  lrn *= pow(clset->steps+1, dec);
    m = (clset->type & CLS_COVARS) ? MAT_UPPER : MAT_DIAG;
    i = clset->clscnt;          /* get the matrix operation mode */
    if (clset->type & CLS_JOINT) {
      while (--i > 0)           /* if to combine (co)variances */
        mat_addx(c[0].cov, c[0].cov, 1, c[i].cov, m);
      i = 1;                    /* combine them in first cluster and */
    }                           /* then process only this cluster */
    for (c += i; --i >= 0; ) {  /* traverse the clusters */
      w = mat_getwgt((--c)->cov);
      if (w >= 1) { eta = lrn/w; dec = 1 -lrn;   }
      else        { eta = lrn;   dec = 1 -lrn*w; }
                                /* compute learning rate and compl. */
      if (clset->type & (CLS_COVARS|CLS_VARS)) {
        mat_mulsx(c->cov, c->cov, eta, m);
        if (dec > 0) mat_addx(c->cov, c->cov, dec, c->smp, m); }
      else {                    /* if adaptable (co)variances */
        w = eta *(mat_trace(c->cov) /clset->incnt);
        if (dec > 0) w += dec *mat_get(c->smp, 0, 0);
        mat_set(c->cov, 0, 0, w);
      }                         /* if only isoptropic variance, */
    }                           /* use specialized procedure */
  }                             /* (update (co)variances/sizes) */

  /* --- compute new weights --- */
  if (clset->type & CLS_WEIGHT){/* if to use cluster weights */
    dec = clset->decay[2];      /* get the decay parameter and */
    lrn = clset->rates[2];      /* compute the next learning rate */
    if ((dec > 0) && (dec < 1)) lrn *= pow(dec, clset->steps);
    else if          (dec < 0)  lrn *= pow(clset->steps+1, dec);
    dec = 1 -lrn; sum = 0;      /* traverse the clusters */
    for (c += i = clset->clscnt; --i >= 0; ) {
      --c; w = lrn *mat_getwgt(c->cov) +dec *mat_getwgt(c->smp);
      if (!(w >= MINWEIGHT)) w = MINWEIGHT;
      sum += mat_setwgt(c->cov, w);      /* adapt the weights and */
    }                           /* sum them for the normalization */
    sum = 1/sum;                /* compute the normalization factor */
    for (c += i = clset->clscnt; --i >= 0; )
      mat_mulwgt((--c)->cov, sum);
  }                             /* normalize the new cluster weights */

  clset->steps++;               /* count the update step */
}  /* _complrn() */

/*----------------------------------------------------------------------
  Error Backpropagation Function
----------------------------------------------------------------------*/

static void _backprop (CLSET *clset)
{                               /* --- error backpropagation */
  int     i;                    /* loop variable */
  CLUSTER *c;                   /* to traverse the clusters */
  double  lrc, lrv, t, v;       /* learning rates, buffers */
  MATRIX  *mat;                 /* exchange buffer */

  assert(clset);                /* check the function argument */
  lrc = -clset->rates[0];       /* get the center   learning rate */
  lrv = -clset->rates[1];       /* and the (co)var. learning rate */
  for (c = clset->cls +(i = clset->clscnt); --i >= 0; ) {
    --c;                        /* traverse the clusters */
    mat = c->smp; c->smp = c->cov; c->cov = mat;
    c->ctr = mat_vector(c->cov);/* exchange the covariance matrix */
    c->sum = mat_vector(c->smp);/* with the aggregation matrix */
    vec_add(c->ctr, clset->incnt, c->sum, lrc, c->ctr);
    if (clset->type & CLS_COVARS) { /* if adaptable covariances */
      mat_mulsx(c->cov, c->cov, lrv, MAT_UPPER);
      mat_addx (c->cov, c->cov, 1, c->inv, MAT_UPPER);
      if (mat_chinvx(c->inv, c->cov) < 0)
        mat_copyx(c->cov, c->smp, MAT_UPPER); }
    else if (clset->type & CLS_VARS) {
      mat_diamuls(c->cov, lrv); /* if adaptable variances */
      mat_diaadd (c->cov, 1, c->inv);
      if (mat_diainv(c->cov, c->cov) < 0)
        mat_copyx(c->cov, c->smp, MAT_DIAG); }
    else if (clset->type & CLS_SIZE) {
      v = mat_get(c->inv, 0,0); /* if adaptable size */