*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 theunnecessary function call for each parameter update as well as theadditional 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 --- */