cluster2.c

聚类算法全集以及内附数据集

    *prv = grd; }               /* note the current gradient */
  else if (t < 0) {             /* if gradients have opposite signs */
    *chg *= clset->shrink;      /* decrease the learning rate */
    if (*chg < clset->minchg) *chg = clset->minchg;
    *prv = 0; }                 /* suppress a change in the next step */
  else {                        /* if one gradient is zero */
    *prv = grd; }               /* only note the current gradient */
  if (grd > 0) return  *chg;    /* return the parameter change */
  if (grd < 0) return -*chg;    /* (use only sign of gradient) */
  return 0;                     /* return 0 if gradient is 0 */
}  /* _resilient() */

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

static double _quick (CLSET *clset,
                      double grd, double *prv, double *chg)
{                               /* --- quickprop analog */
  double t;                     /* temporary buffer */
  double m;                     /* maximal fraction of new derivative */
  
  t = *prv -grd;                /* compute the change of the gradient */
  if      (*chg *t >= 0)        /* if the parabola opens downwards, */
    *chg = grd;                 /* do a standard update step */
  else if (*prv > 0) {          /* if previous gradient was positive */
    m = clset->growth /(clset->growth +1);
    if (grd >= *prv *m)         /* if the growth factor would become */
      *chg *= clset->growth;    /* too large, use the maximal factor */
    else                        /* otherwise compute the factor */
      *chg *= grd /t; }         /* for a jump to the minimum */
  else if (*prv < 0) {          /* if previous gradients was negative */
    m = clset->growth /(clset->growth +1);
    if (grd <= *prv *m)         /* if the growth factor would become */
      *chg *= clset->growth;    /* too large, use the maximal factor */
    else                        /* otherwise compute the factor */
      *chg *= grd /t; }         /* for a jump to the minimum */
  else                          /* if this is the first update */
    *chg = grd;                 /* do a standard update step */
  if      (*chg >  clset->maxchg) *chg =  clset->maxchg; /* clamp */
  else if (*chg < -clset->maxchg) *chg = -clset->maxchg; /* change */
  *prv = grd;                   /* note the new gradient and */
  return *chg;                  /* return the parameter change */
}  /* _quick() */

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

static UPDATEFN *_updatefn[] = {
  /* CLS_NONE       0x00 */  _standard,
  /* CLS_EXPAND     0x10 */  _expand,
  /* CLS_MOMENTUM   0x20 */  _momentum,
  /* CLS_ADPATIVE   0x30 */  _adaptive,
  /* CLS_RESILIENT  0x40 */  _resilient,
  /* CLS_QUICK      0x50 */  _quick,
};                              /* list of update functions */

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

static double _stdupd (CLSET *clset)
{                               /* --- standard update method */
  int     i, n;                 /* loop variables */
  int     type;                 /* cluster type flags */
  CLUSTER *p;                   /* to traverse the clusters */
  double  *s, *c;               /* to access the vectors */
  double  max = 0, t, a = 0;    /* maximal change, buffers */

  assert(clset);                /* check the function arguments */
  type = clset->type;           /* get the cluster type flags */
  for (p = clset->cls +(n = clset->clscnt); --n >= 0; ) {
    --p;                        /* traverse the clusters */
    if (type & CLS_WEIGHT)      /* -- if adaptable weights */
      p->wgt = p->nw;           /* copy the new cluster weight */
    if (p->d2 < 0) continue;    /* skip clusters not to be updated */
    s = p->sum; c = p->ctr;     /* get aggregation vector and center */
    for (i = clset->incnt; --i >= 0; ) {
      t = fabs(s[i] -c[i]); if (t > max) max = t;
      c[i] = s[i];              /* determine the maximal change and */
    }                           /* copy the new center coordinates */
    if (type & CLS_SIZE) {      /* if the cluster size is adaptable, */
      a = 1; p->var = p->msd; } /* copy the new cluster size */
    else {                      /* otherwise (fixed cluster size) */
      a = p->var /p->msd; }     /* compute the rescaling factor */
    p->scl = 1/a;               /* and its reciprocal value */
    if (type & CLS_COVARS) {    /* -- if adaptable covariances */
      mat_chinv(p->inv,p->inv); /* invert the new covariance matrix */
      if (a != 1) {             /* if fixed cluster size */
        mat_trmuls(p->smp, p->smp, MAT_UPPER, a);
        mat_trmuls(p->inv, p->inv, MAT_UPPER, p->scl);
      }                         /* rescale to the old cluster size */
      mat_copy  (p->cov, p->smp, MAT_UPPER|MAT_NOTBUF);
      mat_tr2sym(p->cov, p->cov, MAT_UPPER);  /* copy and complete */
      mat_tr2sym(p->inv, p->inv, MAT_UPPER);  /* the cov. matrices */
      mat_init(p->smp, MAT_ZERO, NULL); }     /* reset aggregation */
    else if (type & CLS_VARS) { /* -- if adaptable variances */
      for (i = clset->incnt; --i >= 0; ) {
        t = a *mat_get(p->smp, i, i);
        if (t < MINVAR) t = MINVAR;
        if (t > MAXVAR) t = MAXVAR;
        mat_set(p->cov, i, i,   t);
        mat_set(p->inv, i, i, 1/t);    /* copy the new variances */
      }                         /* compute the inverse variances */
      mat_init(p->smp, MAT_RESET, NULL);
      mat_diainit(p->smp, -1);} /* reset the aggregation matrix */
    else {                      /* -- if isotropic variance */
      mat_init(p->smp, MAT_RESET, NULL);
      if (type & CLS_SIZE) mat_diainit(p->smp, -1);
    }                           /* reset the aggregation matrix */
  }                             /* and reinitialize the variances */
  return max;                   /* return the maximal change */
}  /* _stdupd() */              /* of a center coordinate */

/*----------------------------------------------------------------------
A special function is used for the standard version to avoid the
unnecessary function call for each parameter update as well as the
additional matrix decompositions needed for a neural network update.
----------------------------------------------------------------------*/

static double _nnupd (CLSET *clset)
{                               /* --- neural network update methods */
  int      i, k, n;             /* loop variables */
  int      type;                /* cluster type flags */
  CLUSTER  *p;                  /* to traverse the clusters */
  double   *s, *c, *x, *b;      /* to access the vectors */
  double   max = 0;             /* maximum change of a parameter */
  double   det, var;            /* determinant and isotropic variance */
  double   t, d, v, u, a;       /* temporary buffers */
  UPDATEFN *update;             /* parameter update function */

  assert(clset);                /* check the function arguments */
  type = clset->type;           /* get the cluster type flags */
  for (p = clset->cls +(n = clset->clscnt); --n >= 0; ) {
    --p;                        /* traverse the clusters */

    /* -- update the cluster weight -- */
    if (type & CLS_WEIGHT)      /* if adaptable cluster weights, */
      p->wgt = p->nw;           /* copy the new cluster weight */
    if (p->d2 < 0) continue;    /* skip clusters not to be updated */

    /* -- update the cluster center -- */
    update = _updatefn[(clset->method & CLS_MODIFIER) >> 4];
    s = p->sum; c = p->ctr;     /* get aggregation vector, center, */
    x = p->chc; b = p->bfc;     /* and the parameter vectors */
    for (i = clset->incnt; --i >= 0; ) {
      d = s[i] -c[i];           /* compute the center gradient */
      t = fabs(d); if (t > max) max = t;
      c[i] += update(clset, d, b+i, x+i);
    }                           /* update the center coordinates */

    /* -- update the (co)variances -- */
    a = (type & CLS_SIZE)       /* compute the rescaling factor */
      ? 1 : p->var /p->msd;     /* (if the cluster size is fixed) */
    p->scl = 1/a;               /* and its reciprocal value */
    if (type & (CLS_COVARS|CLS_VARS)) {
      for (i = clset->incnt; --i >= 0; ) {
        for (k = (type & CLS_COVARS) ? clset->incnt : i+1; --k >= i; ) {
          v = mat_get(p->cov, i, k);    /* get the covariance matrix */
          d = mat_get(p->smp, i, k) *a -v; /* element and its change */
          t = mat_get(p->bfv, i, k);    /* and update its value with */
          u = mat_get(p->chv, i, k);    /* possible add. parameters */
          v += update(clset, d, &t, &u);
          mat_set(p->cov, i, k, v);     /* set the new value of the */
          mat_set(p->bfv, i, k, t);     /* covariance matrix element, */
          mat_set(p->chv, i, k, u);     /* note the updated params. */
        }                       /* (compute a new covariance matrix */
      } }                       /*  or at least new variances) */
    else if (type & CLS_SIZE) { /* -- if isotropic variance */
      t = mat_get(p->bfv, 0,0); /* get the additional parameters */
      u = mat_get(p->chv, 0,0); /* for the variance update */
      p->var += update(clset, p->msd -p->var, &t, &u);
      if (p->var < MINVAR) p->var = MINVAR;
      if (p->var > MAXVAR) p->var = MAXVAR;
      mat_set(p->bfv, 0, 0, t); /* update the isotropic variance, */
      mat_set(p->chv, 0, 0, u); /* ensure that it is positive, */
    }                           /* and note the updated parameters */

    /* -- compute inverses and reinitialize -- */
    if (type & CLS_COVARS) {    /* -- if adaptable covariances */
      if (mat_chdecom(p->inv, p->cov) != 0) {
        mat_copy(p->cov, p->smp, MAT_UPPER);
        mat_chdecom(p->inv, p->cov);
      }                         /* decompose new covariance matrix */
      det = mat_chdet(p->inv);  /* compute the new determinant */
      var = ((det >= MINDET) && (det <= MAXDET))
          ? pow(det,            1.0/clset->incnt)
          : exp(mat_chlogd(p->inv) /clset->incnt);
      if (var < MINVAR) var = MINVAR;
      if (var > MAXVAR) var = MAXVAR;
      mat_chinv(p->inv,p->inv); /* invert the new covariance matrix */
      if (type & CLS_SIZE)      /* if the cluster size is adaptable, */
        p->var = var;           /* copy the new cluster size */
      else {                    /* if the cluster size is fixed */
        mat_trmuls(p->cov, p->cov, MAT_UPPER, p->var /var);
        mat_trmuls(p->inv, p->inv, MAT_UPPER, var /p->var);
      }                         /* scale to the fixed cluster size */
      mat_tr2sym(p->cov, p->cov, MAT_UPPER);  /* complete the      */
      mat_tr2sym(p->inv, p->inv, MAT_UPPER);  /* two matrices and  */
      mat_init(p->smp, MAT_ZERO, NULL); }     /* reset aggregation */
    else if (type & CLS_VARS) { /* -- if adaptable variances */
      det = mat_diaprod(p->cov); /* compute the new determinant */
      var = ((det >= MINDET) && (det <= MAXDET))
          ? pow(det,            1.0/clset->incnt)
          : exp(mat_dialog(p->smp) /clset->incnt);
      if (var < MINVAR) var = MINVAR;
      if (var > MAXVAR) var = MAXVAR;
      if (type & CLS_SIZE) {    /* if the cluster size is adaptable, */
        a = 1; p->var = var; }  /* copy the new cluster size */
      else {                    /* otherwise (fixed cluster size) */
        a = p->var /var; }      /* compute the rescaling factor */
      for (i = clset->incnt; --i >= 0; ) {
        t = fabs(s[i] -c[i]); if (t > max) max = t;
        c[i] = s[i];            /* copy the new center coordinates */
        t = mat_get(p->cov, i, i) *a;
        if (t < MINVAR) t = MINVAR;
        if (t > MAXVAR) t = MAXVAR;
        mat_set(p->cov, i, i, t);
        mat_set(p->inv, i, i, 1/t);    /* copy the variances and */
      }                         /* compute the inverse variances */
      mat_init(p->smp, MAT_RESET, NULL);
      mat_diainit(p->smp, -1);} /* reset the aggregation matrix */
    else {                      /* --- if isotropic variance */
      if (type & CLS_SIZE) {    /* if the cluster size is adaptable, */
        if (p->var < MINVAR) p->var = MINVAR;
        if (p->var > MAXVAR) p->var = MAXVAR;
        mat_diainit(p->smp,-1); /* reinitialize the variances */
      }
      mat_init(p->smp, MAT_RESET, NULL);
    }                           /* reset the aggregation matrix */
  }
  return max;                   /* return the maximal change */
}  /* _nnupd() */               /* of a center coordinate */

/*----------------------------------------------------------------------
  Main Functions
----------------------------------------------------------------------*/

void cls_init (CLSET *clset, int mode, double range,
               double randfn(void), const double *vec)
{                               /* --- initialize a set of clusters */
  int     i, k, n;              /* loop variables */
  CLUSTER *p;                   /* to traverse the clusters */
  NSTATS  *nst;                 /* simple numerical statistics */
  double  *c, *b;               /* cluster center, buffer */
  double  w, x, d, m;           /* cluster weight, buffers */

  assert(clset && randfn);      /* check the function arguments */

  /* --- initialize other parameters --- */
  if (((mode & 0xf) != CLS_POINTS)    /* if not called second time */
  ||  (clset->init <= 0) || (clset->init >= clset->clscnt)) {
    w = 1.0 /clset->clscnt;     /* initial relative cluster weight */
    for (p = clset->cls +(k = clset->clscnt); --k >= 0; ) {
      (--p)->wgt = w;           /* traverse the clusters and */
      p->var = p->scl = p->d2 = 1;   /* init. the parameters */
      mat_init(p->cov, MAT_UNIT, NULL);
      mat_init(p->inv, MAT_UNIT, NULL);
      mat_init(p->smp, MAT_ZERO, NULL);
    }                           /* initialize the matrices */
  }                             /* for all clusters */

  /* --- special initialization for centers at origin --- */
  if (mode & CLS_ORIGIN) {      /* if to set centers to origin */
    for (p = clset->cls +(k = clset->clscnt); --k >= 0; ) {
      for (c = (--p)->ctr +(i = clset->incnt); --i >= 0; ) {
        *--c = 0;               /* set all coordinates to zero */
        mat_set(p->cov, i, i, x = 0.4 +0.8*randfn());
        mat_set(p->inv, i, i, 1/x);
      }                         /* initialize the diagonal of */
    } return;                   /* the covariance matrix randomly */
  }                             /* and set the inverse matrix */

  /* --- initialize cluster centers --- */