cluster1.c

it is the Data Mining Algorithm source code.

The above functions do a safe decomposition of the covariance matrix or
a safe computation of the isotropic variance of a cluster by shifting
the eigenvalues in case of a failure. This is repeated with larger
and larger shifts and only if after a shift by a fairly large value
(MINVAR *2^96) the decomposition/variance computation still fails, a
unit matrix is set for the covariance matrix (and its decomposition).
----------------------------------------------------------------------*/

void cls_vars (CLSET *clset, int eigen)
{                               /* --- compute isotropic variances */
  int     i;                    /* loop variable */
  CLUSTER *c;                   /* to traverse the clusters */

  assert(clset);                /* check the function argument */
  clset->eigen = 0;             /* clear the eigenvalue decomp. flag */
  i = (clset->type & CLS_JOINT) ? 1 : clset->clscnt;
  c = clset->cls +i;            /* get the cluster array */
  if (clset->type & CLS_COVARS){/* if full covariance matrix */
    if ((clset->fwexp >  0)     /* if feature weighting */
    &&  (clset->fwexp != 1))    /* (exponent or alternative), */
      eigen = -1;               /* eigenvalue decomposition is needed */
    clset->eigen = eigen;       /* note the eigenvalue decomp. flag */
    while (--i >= 0)            /* traverse the clusters and */
      _covmat(--c, eigen); }    /* decompose the covariance matrices */
  else if (clset->type & CLS_VARS) {   /* if only variances, */
    while (--i >= 0)            /* traverse the clusters and */
      _diamat(--c); }           /* handle diagonal matrices */
  else {                        /* if only isotropic variances */
    while (--i >= 0) { --c;     /* traverse the clusters */
      c->var = mat_get(c->cov, 0, 0);
      if (!(c->var >= MINVAR)) mat_set(c->cov, 0, 0, c->var = MINVAR);
      if (!(c->var <= MAXVAR)) mat_set(c->cov, 0, 0, c->var = MAXVAR);
    }                           /* get the isotropic variance and */
  }                             /* clamp it to a reasonable range */
}  /* cls_vars() */

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

void cls_gauge (CLSET *clset)
{                               /* --- measure the cluster sizes */
  int     i;                    /* loop variable */
  CLUSTER *c;                   /* to traverse the clusters */
  double  t;                    /* buffer for computations */

  assert(clset);                /* check the function argument */
  i = (clset->type & CLS_JOINT) ? 1 : clset->clscnt;
  c = clset->cls +i;            /* get the cluster array */
  if ((clset->fwexp <= 0)       /* if no feature weighting, */
  ||  (clset->fwexp == 1)       /* the isotropic variance is the size */
  || !(clset->type & (CLS_COVARS|CLS_VARS))) {
    while (--i >= 0) { --c;     /* traverse the clusters and */
      c->scale = c->var; } }    /* copy the isotropic variance */
  else {                        /* if feature weighting */
    assert(clset->eigen         /* check for eigenvalue decomposition */
      || !(clset->type & CLS_COVARS));
    while (--i >= 0) { --c;     /* traverse the clusters */
      t = (clset->type & CLS_COVARS)
        ? vec_sumrec(c->dif, clset->incnt) : mat_diarec(c->cov);
      c->scale = (t > MINVAR) ? 1/t : MAXVAR;
    }                           /* compute the cluster size */
  }                             /* the sum of the feature weights */
}  /* cls_gauge() */

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

void cls_resize (CLSET *clset, int decomp)
{                               /* --- (re)scale the clusters */
  int     i;                    /* loop variable */
  CLUSTER *c;                   /* to traverse the clusters */

  assert(clset);                /* check the function argument */
  i = (clset->type & CLS_JOINT) ? 1 : clset->clscnt;
  c = clset->cls +i;            /* get the cluster array */
  if (clset->type & CLS_SIZE) { /* if adaptable cluster size, */
    clset->resized = 1;         /* set flag for resized clusters */
    while (--i >= 0) { --c;     /* traverse the clusters */
      c->size  = c->scale;      /* simple copy the cluster size */
      c->scale = 1;             /* and set a unit scaling factor */
    } }
  else {                        /* if the cluster size is fixed */
    while (--i >= 0) { --c;     /* traverse the clusters */
      c->scale = c->size /c->scale;  /* compute the scaling factor */
      if ((c->scale == 1) || !(clset->type & CLS_COVARS) || !decomp)
        continue;               /* check whether scaling is needed */
      if (clset->eigen) vec_muls(c->dif,clset->incnt, c->dif, c->scale);
      else mat_mulsx(c->inv, c->inv, sqrt(c->scale), MAT_LOWER);
    }                           /* rescale the decomposition */
  }                             /* of the covariance matrix */
  if (!(clset->type & (CLS_COVARS|CLS_VARS)) && clset->resized) {
    i = (clset->type & CLS_JOINT) ? 1 : clset->clscnt;
    for (c += i; --i >= 0; ) {  /* traverse the clusters */
      --c; mat_diasetx(c->cov, mat_get(c->cov, 0, 0)); }
  }                             /* (re)set the diagonal */
}  /* cls_resize() */

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

void cls_ftwgts (CLSET *clset)
{                               /* --- compute feature weights */
  int     i, k, n;              /* loop variables, number of dims. */
  CLUSTER *c;                   /* to traverse the clusters */
  double  x, v, sum;            /* buffers for computations */
  double  off;                  /* transformation offset */

  assert(clset);                /* check the function argument */
  i = (clset->type & CLS_JOINT) ? 1 : clset->clscnt;
  c = clset->cls +i;            /* get the cluster array and */
  n = clset->incnt;             /* the number of dimensions */
  x = clset->fwexp;             /* get the feature weight exponent */

  /* --- compute feature weights --- */
  if ((x > 0) && (x != 1)       /* if to do feature weighting */
  &&  (clset->type & (CLS_COVARS|CLS_VARS))) {
    assert(clset->eigen         /* check for eigenvalue decomposition */
      || !(clset->type & CLS_COVARS));
    while (--i >= 0) { --c;     /* traverse the clusters */
      sum = 1/c->scale;         /* compute the scaling factor */
      for (k = n; --k >= 0; ) { /* traverse the dimensions */
        v = (clset->type & CLS_COVARS)
          ? c->dif[k] : mat_get(c->cov, k, k);
        c->dif[k] = sum / ((v > 0) ? v : MINVAR);
      }                         /* compute the feature weights */
      if (x > 1) continue;      /* if normal exponent, skip extension */
      vec_copy (clset->buf, c->dif, n);
      v_dblsort(clset->buf, n); /* copy the inverse eigenvalues */
      v_dblrev (clset->buf, n); /* and sort them descendingly */
      for (sum = v = clset->buf[k = 0]; ++k < n; ) {
        v = clset->buf[k];      /* traverse the sorted eigenvalues */
        if (v *(1 +x*k) <= x *(sum +v)) break;
        sum += v;               /* find the number of features */
      }                         /* that receive a positive weight */
      /* Note that if two (or more) degrees of membership are equal, */
      /* then either all of them or none of them satisfy the above   */
      /* condition for a non-zero transformed degree of membership.  */
      if (!(sum > 0)) sum = 1;  /* make sure the sum does not vanish */
      off = x *(v = 1/(1-x));   /* compute transformation parameters */
      sum = (1 +x *(k-1)) *v /sum;
      for (k = n; --k >= 0; ) { /* traverse the inverse eigenvalues */
        v = sum *c->dif[k] -off;
        c->dif[k] = (v > 0) ? v : 0;
      }                         /* compute the feature weights, */
    } }                         /* possibly setting some to zero */

  /* --- compute inverse eigenvalues --- */
  else if (clset->eigen) {      /* if eigenvalue decomposition */
    while (--i >= 0) { --c;     /* traverse the clusters */
      sum = 1/c->scale;         /* compute the scaling factor */
      for (k = n; --k >= 0; ) {
        v = c->dif[k];          /* traverse the eigenvalues */
        c->dif[k] = sum / ((v > 0) ? v : MINVAR);
      }                         /* compute inverse eigenvalues */
    } }                         /* in a safe way */
  else if (clset->type & CLS_COVARS)
    { }                         /* if Cholesky decomp., do nothing */
  else if (clset->type & CLS_VARS) {
    while (--i >= 0) { --c;     /* traverse the clusters */
      sum = 1/c->scale;         /* compute the scaling factor */
      for (k = n; --k >= 0; ) { /* traverse the variances */
        v = mat_get(c->cov, k, k);
        c->dif[k] = sum / ((v > 0) ? v : MINVAR);
      }                         /* compute inverse variances */
    } }                         /* in a safe way */
  else if ((clset->type & CLS_SIZE) || clset->resized) {
    while (--i >= 0) { --c;     /* traverse the clusters */
      v = mat_get(c->cov, 0, 0);
      c->dif[0] = (1/c->scale) /((v > 0) ? v : MINVAR);
    }                           /* invert the isotropic variance */ 
  }                             /* in a safe way */
}  /* cls_ftwgts() */

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

void cls_invert (CLSET *clset)
{                               /* --- compute inverse matrices */
  int     i, k, n;              /* loop variables, buffers */
  CLUSTER *c;                   /* to traverse the clusters */
  double  t, x, b;              /* buffers for computations */

  assert(clset);                /* check the function argument */
  i = (clset->type & CLS_JOINT) ? 1 : clset->clscnt;
  c = clset->cls +i;            /* get the cluster array */
  x = clset->fwexp;             /* get the feature weight exponent */

  /* --- if feature weights --- */
  if ((x > 0) && (x != 1)       /* if to do feature weighting */
  &&  (clset->type & (CLS_COVARS|CLS_VARS))) {
    assert(clset->eigen         /* check for eigenvalue decomposition */
      || !(clset->type & CLS_COVARS));
    b = (x < 1) ? (1-x)/(1+x) : 0; /* compute transform. parameter */
    n = clset->incnt;           /* get the number of dimensions */
    while (--i >= 0) { --c;     /* traverse the clusters */
      k = n;                    /* if normal feature weight exponent, */
      if (x > 1) {              /* simply adapt inverse eigenvalues */
        if (x == 2) while (--k >= 0) c->dif[k] *= c->dif[k];
        else        while (--k >= 0) c->dif[k]  = pow(c->dif[k], x); }
      else {                    /* if alternative to an exponent */
        while (--k >= 0) {      /* traverse the dimensions */
          t = (b *(c->dif[k] -1) +1) *c->dif[k];
          c->dif[k] = (t > 0) ? t : 0;
	}                       /* transform the feature weights */
      }                         /* with alternative formula */
      if (!(clset->type & CLS_COVARS)) mat_diaset(c->inv, c->dif);
      else mat_mulmdmt(c->inv, c->inv, c->dif);
    } }                         /* compute the distance norm matrix */

  /* --- if normal (co)variances --- */
  else {                        /* if no feature weighting */
    if (clset->eigen)           /* if eigenvalue decomposition */
      while (--i >= 0) { --c;   /* multiply eigenvalues and -vectors */
        mat_mulmdmt(c->inv, c->inv, c->dif); }
    else if (clset->type & CLS_COVARS) {
      while (--i >= 0) { --c;   /* if Cholesky decomposition */
        mat_tr2sym(c->cov, c->cov, MAT_UPPER);
        mat_chinv (c->inv, c->inv);
        mat_tr2sym(c->inv, c->inv, MAT_UPPER);
      } }                       /* compute inverse covariance matrix */
    else if (clset->type & CLS_VARS)
      while (--i >= 0) { --c;   /* if only variances, set inverses */
        mat_diaset(c->inv, c->dif); }
    else if (clset->resized) {  /* if clusters were resized, */
      while (--i >= 0) { --c;   /* set the inverse variance */
        mat_diasetx(c->inv, c->dif[0]); }
    }                           /* (which is also stored only in */
  }                             /* c->dif first for consistency) */
  clset->resized = 0;           /* clear flag for resized clusters */
}  /* cls_invert() */

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

void cls_reinit (CLSET *clset)
{                               /* --- reinit. aggregation matrices */
  int     i, m;                 /* loop variable, buffer */
  CLUSTER *c;                   /* to traverse the clusters */

  if      (clset->type & CLS_COVARS)          m = MAT_UPPER;
  else if (clset->type & (CLS_VARS|CLS_SIZE)) m = MAT_DIAG;
  else                                        m = 0;
  m |= MAT_VECTOR|MAT_WEIGHT;   /* get the initialization mode */
  for (c = clset->cls +(i = clset->clscnt); --i >= 0; ) {
    (--c)->dec = 0;             /* traverse the clusters */
    mat_init(c->smp, m, NULL);  /* clear the decomposition flag */
  }                             /* and reinitialize the matrices */
}  /* cls_reinit() */

/*----------------------------------------------------------------------
  Cluster Set Functions
----------------------------------------------------------------------*/

CLSET* cls_create (int incnt, int clscnt)
{                               /* --- create a set of clusters */
  int     i;                    /* loop variable */
  CLSET   *clset;               /* created set of clusters */
  CLUSTER *c;                   /* to traverse the clusters */

  assert((incnt  > 0)           /* check the function arguments */