/*----------------------------------------------------------------------  File    : cluster2.c  Contents: cluster and cluster set management (update functions)  Author  : Christian Borgelt  History : 05.09.2001 file created as cluster1.c            12.09.2001 function cls_init completed            15.09.2001 first version of fuzzy c-means completed            16.09.2001 hard c-means algorithm added (msexp <= 0)            09.09.2002 neural network update methods added            29.01.2003 some cleanup of neural network methods            31.01.2003 initialization for resilient method changed            18.02.2003 multiple init. in mode CLS_POINTS removed            07.06.2003 cluster size adaptation added            29.10.2003 bug in function cls_init fixed (CLS_POINTS)            22.02.2004 size computation for spherical clusters modified            23.02.2004 bug in function _quick fixed            29.02.2004 some bugs in cluster size computation fixed            01.03.2004 update functions moved to this file            02.03.2004 shape and size regularization added            19.03.2004 weight/prior regularization added            12.04.2004 competitive learning function completed            13.04.2004 regularization adapted for competitive learning            14.04.2004 treatment of (almost) empty clusters improved            15.04.2004 update function improved (loop moved)            23.04.2004 some numeric problems with high-dim. data solved            26.04.2004 more numeric problems with high-dim. data solved            27.04.2004 upper learning rate bound added (alt. est. step)            30.04.2004 rescaling removed from parameter est. functions            13.07.2004 normalization of center vectors added            15.07.2004 bug in center normalization for _complrn fixed            28.07.2004 update of centers only added to cls_update            14.08.2004 bug in initialization with CLS_POINTS fixed            18.08.2004 first version of backpropagation functions added----------------------------------------------------------------------*/#include <stdlib.h>#include <float.h>#include <math.h>#include <assert.h>#include "cluster.h"/*----------------------------------------------------------------------  Preprocessor Definitions----------------------------------------------------------------------*/#define MINVAR       1e-12      /* minimal variance */#define MAXVAR       1e+12      /* maximal variance */#define MINDET       1e-48      /* minimal determinant */#define MAXDET       1e+48      /* maximal determinant */#define MINWEIGHT    1e-6       /* 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----------------------------------------------------------------------*/static double _decom (CLUSTER *p){                               /* --- decompose covariance matrix */  int    i;                     /* loop variables */  double t;                     /* buffer fro shift value */  for (t = MINVAR, i = 40; --i >= 0; t += t) {    if (mat_chdecom(p->inv, p->smp) == 0) break;    mat_diaadds(p->smp, t);     /* decompose the covariance matrix */  }                             /* and on failure shift eigenvalues */  if (i < 0) {                  /* if decomposition fails totally */    mat_init(p->smp, MAT_UNIT|MAT_NOTBUF, NULL);    mat_init(p->inv, MAT_UNIT|MAT_NOTBUF, NULL);  }                             /* set a unit covariance matrix */  return mat_chdet(p->inv);     /* return the determinant */}  /* _decom() *//*--------------------------------------------------------------------*/static void _normctr (CLSET *clset, int sum){                               /* --- normalize center vectors */  int     i, k;                 /* loop variables */  CLUSTER *p;                   /* to traverse the clusters */  double  *c;                   /* to traverse the center vector */  double  len;                  /* length of a center vector */    assert(clset);                /* check the function argument */  for (p = clset->cls +(i = clset->clscnt); --i >= 0; ) {    if ((--p)->d2 < 0) continue;/* traverse marked clusters */    c = (sum) ? p->sum : p->ctr;/* get the vector to normalize */    len = 0;                    /* initialize the length */    for (c += k = clset->incnt; --k >= 0; ) {      --c; len += *c * *c; }    /* sum the squared coordinates */    len = sqrt(len);            /* compute the vector length */    len = (len > 0)             /* and the normalization factor */        ? (sum ? mat_weight(p->smp) : 1.0) /len : 1;    if (fabs(1.0 -len) < 1e-12) continue;    for (c += k = clset->incnt; --k >= 0; )      *--c *= len;              /* normalize the center vectors */  }                             /* to unit length */}  /* _normctr() *//*--------------------------------------------------------------------*/static void _zeroctr (CLSET *clset, int sum){                               /* --- set center vectors to origin */  int     i, k;                 /* loop variables */  CLUSTER *p;                   /* to traverse the clusters */  double  *c;                   /* to traverse the center vector */  assert(clset);                /* check the function argument */  for (p = clset->cls +(i = clset->clscnt); --i >= 0; ) {    if ((--p)->d2 < 0) continue;/* traverse marked clusters */    c = (sum) ? p->sum : p->ctr;/* get the vector to zero */    for (c += k = clset->incnt; --k >= 0; )      *--c = 0;                 /* set all coordinates to zero */  }                             /* (move vector to origin) */}  /* _zeroctr() *//*----------------------------------------------------------------------  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     n;                    /* loop variable */  int     type;                 /* cluster type flags */  CLUSTER *p;                   /* to traverse the clusters */  double  det;                  /* determinant of covariance matrix */  assert(clset);                /* check the function argument */  if (clset->method & CLS_ORIGIN)  /* if cluster centers at origin, */    _zeroctr(clset, 1);            /* zero the center vectors */  if (clset->method & CLS_UNIT)    /* if centers on unit sphere, */    _normctr(clset, 1);            /* normalize the center vectors */  type = clset->type;           /* get the cluster type flags */  for (p = clset->cls +(n = clset->clscnt); --n >= 0; ) {    (--p)->nw *= clset->msd[1]; /* normalize cluster weights to sum 1 */    if (p->d2 < 0) continue;    /* skip clusters not to be updated */    if (type & CLS_COVARS) {    /* -- if adaptable covariances */      mat_covar(p->smp, p->smp, 1); /* compute new covariances */      det    = _decom(p);           /* and decompose the matrix */      p->msd = ((det >= MINDET) && (det <= MAXDET))             ? pow(det,            1.0/clset->incnt)             : exp(mat_chlogd(p->inv) /clset->incnt); }    else if (type & CLS_VARS) { /* -- if adaptable variances */      mat_var(p->smp, p->smp, 1);   /* compute new variances */      det    = mat_diaprod(p->smp); /* and the new determinant */      p->msd = ((det >= MINDET) && (det <= MAXDET))             ? pow(det,            1.0/clset->incnt)             : exp(mat_dialog(p->smp) /clset->incnt); }    else if (type & CLS_SIZE)   /* -- if adaptable isotropic variance */      p->msd = mat_isovar(p->sum, p->smp, 1);    else {                      /* -- if fixed isotropic variance */      mat_mean(p->sum, p->smp); /* compute only the mean values */      p->msd = p->var;          /* (i.e. a new cluster center) */    }                           /* copy the old isotropic variance */    if      (p->msd < MINVAR) p->msd = MINVAR;    else if (p->msd > MAXVAR) p->msd = MAXVAR;  }                             /* clamp the variance */}  /* _altopt() *//*----------------------------------------------------------------------  Competitive Learning Function----------------------------------------------------------------------*/static void _complrn (CLSET *clset){                               /* --- competitive learning */  int     i, n;                 /* loop variables */  int     type;                 /* cluster type flags */  CLUSTER *p;                   /* to traverse the clusters */  double  *s, *c;               /* to access the vectors */  double  lrc, lrv = 0, lrw = 0;/* learning rates */  double  eta, dec;             /* learning rate derivates */  double  t, d;                 /* temporary buffers */  assert(clset);                /* check the function argument */  /* --- compute learning rates --- */  t   = clset->decays[0];       /* get the decay parameter and */  lrc = clset->lrates[0];       /* compute the next learning rate */  if ((t > 0) && (t < 1)) lrc *= pow(t, clset->steps);  else if        (t < 0)  lrc *= pow(clset->steps+1, t);  type = clset->type;           /* get the cluster type flags */  if (type & (CLS_COVARS|CLS_VARS|CLS_SIZE)) {    t   = clset->decays[1];     /* get the decay parameter and */    lrv = clset->lrates[1];     /* compute the next learning rate */    if ((t > 0) && (t < 1)) lrv *= pow(t, clset->steps);    else if        (t < 0)  lrv *= pow(clset->steps+1, t);  }                             /* (learning rate for (co)variances) */  if (type & CLS_WEIGHT) {      /* if adaptable weights */    t   = clset->decays[2];     /* get the decay parameter and */    lrw = clset->lrates[2];     /* compute the next learning rate */    if ((t > 0) && (t < 1)) lrw *= pow(t, clset->steps);    else if        (t < 0)  lrw *= pow(clset->steps+1, t);  }                             /* (learning rate for weights) */  clset->steps++;               /* count the update step */  /* --- compute new parameters --- */  for (p = clset->cls +(n = clset->clscnt); --n >= 0; ) {    (--p)->nw *= clset->msd[1]; /* normalize cluster weights to sum 1 */    if (type & CLS_WEIGHT)      /* and compute new cluster weights */      p->nw = (1 -lrw) *p->wgt +lrw *p->nw;    if (p->d2 < 0) continue;    /* skip clusters not to be updated */    t = mat_weight(p->smp);     /* get the aggregation weight and */    eta = (t < 1) ? lrc : lrc/t;/* compute the learning rate and */    s = p->sum; c = p->ctr;     /* get aggregation vector and center */    for (i = clset->incnt; --i >= 0; )      s[i] = c[i] +eta *s[i];   /* compute new center coordinates */    /* Here p->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 p->ctr.         */    if (t < 1) { eta = lrv;   dec = 1 -lrv*t; }    else       { eta = lrv/t; dec = 1 -lrv;   }    dec *= p->scl;              /* compute learning rate and compl. */    if (type & CLS_COVARS) {    /* -- if adaptable covariances */      mat_trmuls(p->smp, p->smp, MAT_UPPER, eta);      if (dec > 0) mat_addx(p->smp, p->smp, dec, p->cov, MAT_UPPER);      d      = _decom(p);       /* update the covariance matrix */      p->msd = ((d >= MINDET) && (d <= MAXDET))             ? pow(d,              1.0/clset->incnt)             : exp(mat_chlogd(p->inv) /clset->incnt); }    else if (type & CLS_VARS) { /* -- if adaptable variances */      for (d = 1, i = clset->incnt; --i >= 0; ) {        t = ((dec > 0) ? dec *mat_get(p->cov, i, i) : 0)                       + eta *mat_get(p->smp, i, i);        if      (t < MINVAR) t = MINVAR;        else if (t > MAXVAR) t = MAXVAR;        d *= mat_set(p->smp, i, i, t);      }                         /* update the covariance matrix */      p->msd = ((d >= MINDET) && (d <= MAXDET))             ? pow(d,              1.0/clset->incnt)             : exp(mat_dialog(p->smp) /clset->incnt); }    else if (type & CLS_SIZE) { /* -- if adaptable isotropic var. */      d = mat_diaprod(p->smp);  /* update the isotropic variance */      t = ((d >= MINDET) && (d <= MAXDET))        ? pow(d,              1.0/clset->incnt)        : exp(mat_dialog(p->smp) /clset->incnt);      p->msd = ((dec > 0) ? dec *p->var : 0) +eta *t; }    else p->msd = p->var;       /* copy the isotropic variance */    if      (p->msd < MINVAR) p->msd = MINVAR;    else if (p->msd > MAXVAR) p->msd = MAXVAR;  }                             /* clamp the variance */  if (clset->method & CLS_ORIGIN)  /* if cluster centers at origin, */    _zeroctr(clset, 1);            /* zero the center vectors */  if (clset->method & CLS_UNIT)    /* if centers on unit sphere, */    _normctr(clset, 1);            /* normalize the center vectors */}  /* _complrn() *//*----------------------------------------------------------------------  Error Backpropagation Function----------------------------------------------------------------------*/static void _backprop (CLSET *clset){                               /* --- error backpropagation */  int     i, n;                 /* loop variables */  int     type;                 /* cluster type flags */  CLUSTER *p;                   /* to traverse the clusters */  double  *s, *c;               /* to access the vectors */  double  lrc, lrv, d, t;       /* learning rates, buffers */