regress.c

it is regression Algorithm in C/C++.

#endif
/*--------------------------------------------------------------------*/

int reg_desc (REGRESS *reg, FILE *file, int mode, int maxlen)
{                               /* --- describe a regression function */
  int  i, k;                    /* loop variables */
  int  *e      = NULL;          /* table of exponents */
  char *indent = "";            /* indentation string */
  #ifdef REG_EXTFN
  char buf[AS_MAXLEN+1];        /* output buffer for target name */
  #endif

  /* --- print header (as a comment) --- */
  if (mode & REG_TITLE) {       /* if the title flag is set */
    i = k = (maxlen > 0) ? maxlen -2 : 70;
    fputs("/*", file); while (--i >= 0) putc('-', file);
    fprintf(file, "\n  regression polynomial\n");
    while (--k >= 0) putc('-', file); fputs("*/\n", file);
  }                             /* print a title header */
  if (maxlen <= 0) maxlen = INT_MAX;

  /* --- print polynomial --- */
  #ifdef REG_EXTFN              /* if to compile extended functions */
  if (reg->attmap) {            /* if based on an attribute set, */
    fputs("polynomial(", file); /* formally enclose the description */
    sc_format(buf, att_name(reg->trgatt), 0);
    fputs(buf, file);           /* print target attribute name */
    fputs(") = {\n", file);     /* and start the description */
    indent = "  ";              /* indent by two characters */
  } else
  #endif                        /* print the basic parameters */
  fprintf(file, "%svars   = %d;\n", indent, reg->dim);
  fprintf(file, "%sdegree = %d;\n", indent, reg->deg);
  if (reg->max > 0)             /* print maximum for logit transform. */
    fprintf(file, "%slogit  = %g;\n", indent, reg->max);
  if (mode & REG_EXPOS) {       /* if to print exponents */
    e = reg->tab +(i = reg->dim-1);
    *--e = reg->deg +1;         /* init. the table of exponents to */
    while (--i > 0) *--e = 0;   /* x_n^(deg+1), i.e. the "previous" */
  }                             /* state from which to count down */
  fprintf(file, "%scoeffs = {", indent);
  for (i = 0; i < reg->cnt; i++) { /* list the coefficients */
    if (i > 0)         fputc(',', file);
    if (!e && (i % 3)) fputc(' ', file);
    else fprintf(file, "\n%s  ", indent);
    fprintf(file, "%+.16g", reg->cfs[i]);
    if (!e) continue;           /* check whether to print exponents */
    if (e[0] > 0) {             /* update the table of exponents: */
      e[1]++; e[0]--; }         /* if first exponent is not zero, */
    else {                      /* shift a 1 from first to second */
      for (k = 0; e[++k] <= 0; );   /* find first nonzero expo. */
      e[0] = e[k]-1; e[k] = 0; e[k+1]++;
    }                           /* move exponents back to first */
    fputs(" /*", file);         /* start the list of exponents */
    for (k = 0; ++k < reg->dim; ) fprintf(file, " %d", e[k-1]);
    fputs(" */", file);         /* terminate the list of exponents */
  }
  fputs(" };\n", file);         /* terminate the coefficient list */
  #ifdef REG_EXTFN              /* if to compile extended functions */
  if (reg->attmap) fputs("};\n", file);
  #endif                        /* terminate the description */
  return ferror(file) ? -1 : 0; /* return the write status */
}  /* reg_desc() */

/*--------------------------------------------------------------------*/
#ifdef REG_PARSE

static int _getint (SCAN *scan, char *key1, char *key2, int min)
{                               /* --- parse an integer number */
  int  num;                     /* parsed number */
  char *s;                      /* end pointer for conversion */

  assert(scan);                 /* check the function argument */
  if ((sc_token(scan) != T_ID)  /* check for an identifier */
  || ((strcmp(sc_value(scan), key1) != 0)
  &&  (strcmp(sc_value(scan), key2) != 0)))
    ERR_STR(key1);              /* check for keyword */
  GET_TOK();                    /* consume keyword */
  GET_CHR('=');                 /* consume '=' */
  if (sc_token(scan) != T_NUM) ERROR(E_NUMEXP);
  num = (int)strtol(sc_value(scan), &s, 0);
  if (*s || (num < min))       ERROR(E_NUMBER);
  GET_TOK();                    /* get and consume the number */
  GET_CHR(';');                 /* consume ';' */
  return num;                   /* return the parsed number */
}  /* _getint() */

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

static int _parse (REGRESS *reg, SCAN *scan)
{                               /* --- parse a regression polynomial */
  int  i;                       /* loop variable */
  char *s;                      /* end pointer for conversion */

  assert(reg && scan);          /* check the function arguments */
  if ((sc_token(scan) == T_ID)  /* check for 'logit' */
  &&  ((strcmp(sc_value(scan), "logit") == 0)
  ||   (strcmp(sc_value(scan), "max")   == 0))) {
    GET_TOK();                  /* consume 'logit' */
    GET_CHR('=');               /* consume '=' */
    if (sc_token(scan) != T_NUM) ERROR(E_NUMEXP);
    reg->max = strtod(sc_value(scan), &s);
    if (*s || (reg->max <= 0))   ERROR(E_NUMBER);
    GET_TOK();                  /* get and consume the num. of vars. */
    GET_CHR(';');               /* consume ';' */
  }
  if ((sc_token(scan) != T_ID)  /* check for an identifier */
  || ((strcmp(sc_value(scan), "coefficients") != 0)
  &&  (strcmp(sc_value(scan), "coeffs")       != 0)))
    ERR_STR("coeffs");          /* check for 'coeffs' */
  GET_TOK();                    /* consume 'coeffs' */
  GET_CHR('=');                 /* consume '=' */
  GET_CHR('{');                 /* consume '{' */
  for (i = 0; i < reg->cnt; i++) {
    if (i > 0) {GET_CHR(',');}  /* traverse the coefficients */
    if (sc_token(scan) != T_NUM) ERROR(E_NUMEXP);
    reg->cfs[i] = strtod(sc_value(scan), &s);
    if (*s)                      ERROR(E_NUMBER);
    GET_TOK();                  /* get and consume */
  }                             /* a regression coefficient */
  GET_CHR('}');                 /* consume '}' */
  GET_CHR(';');                 /* consume ';' */
  return 0;                     /* return 'ok' */
}  /* _parse() */

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

REGRESS* reg_parse (SCAN *scan)
{                               /* --- parse a regression polynomial */
  int     dim;                  /* number of dimensions/variables */
  int     deg;                  /* degree of the polynomial */
  REGRESS *reg;                 /* parsed regression polynomial */

  assert(scan);                 /* check the function argument */
  pa_init(scan);                /* initialize parsing */
  dim = _getint(scan, "vars",   "dim", 2);
  if (dim < 0) return NULL;     /* get the number of variables */
  deg = _getint(scan, "degree", "deg", 1);
  if (deg < 0) return NULL;     /* get the degree of the polynomial */
  reg = reg_create(dim, deg);   /* create a regression object */ 
  if (!reg)    return NULL;
  if (_parse(reg, scan) != 0) { /* parse the regression polynomial */
    reg_delete(reg); return NULL; }
  return reg;                   /* return the parsed polynomial */
}  /* reg_parse() */

/*--------------------------------------------------------------------*/
#ifdef REG_EXTFN

static int _header (SCAN *scan, ATTSET *attset)
{                               /* --- parse header */
  int t;                        /* token, target identifier */

  assert(attset && scan);       /* check the function arguments */
  if ((sc_token(scan) != T_ID)
  || ((strcmp(sc_value(scan), "polynomial") != 0)
  &&  (strcmp(sc_value(scan), "poly")       != 0)
  &&  (strcmp(sc_value(scan), "regression") != 0)
  &&  (strcmp(sc_value(scan), "regress")    != 0)
  &&  (strcmp(sc_value(scan), "reg")        != 0)))
    ERR_STR("polynomial");      /* check for 'polynomial' */
  GET_TOK();                    /* consume 'polynomial' */
  GET_CHR('(');                 /* consume '(' */
  t = sc_token(scan);           /* check for an attribute name */
  if ((t != T_ID) && (t != T_NUM)) ERROR(E_ATTEXP);
  t = as_attid(attset, sc_value(scan));
  if (t < 0)                       ERROR(E_UNKATT);
  GET_TOK();                    /* consume the attribute name */
  GET_CHR(')');                 /* consume ')' */
  GET_CHR('=');                 /* consume '=' */
  GET_CHR('{');                 /* consume '{' */
  return t;                     /* return target identifier */
}  /* _header() */

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

static int _trailer (SCAN *scan)
{                               /* --- parse trailer */
  assert(scan);                 /* check the function argument */
  GET_CHR('}');                 /* consume '}' */
  GET_CHR(';');                 /* consume ';' */
  return 0;                     /* return 'ok' */
}  /* _trailer() */

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

REGRESS* reg_parsex (SCAN *scan, ATTMAP *attmap)
{                               /* --- parse a regression polynomial */
  int     trgid;                /* identifier of target attribute */
  int     deg;                  /* degree of the polynomial */
  REGRESS *reg;                 /* parsed multivariate polynomial */

  assert(scan && attmap);       /* check the function arguments */
  pa_init(scan);                /* initialize parsing */
  trgid = _header(scan, am_attset(attmap));
  if (trgid < 0) return NULL;   /* parse the header and */
  am_target(attmap, trgid);     /* set the target attribute */
  deg = _getint(scan, "degree", "deg", 1);
  if (deg   < 0) return NULL;   /* get the degree of the polynomial */
  reg = reg_createx(attmap,deg);/* and create a regression object */ 
  if (!reg) { sc_error(scan, E_NOMEM); return NULL; }
  if ((_parse(reg, scan) != 0)  /* parse the polynomial */
  ||  (_trailer(scan)    != 0)) {    /* and the trailer */
    reg_delete(reg); return NULL; }
  return reg;                   /* return the parsed polynomial */
}  /* reg_parsex() */

#endif  /* #ifdef REG_EXTFN */
#endif  /* #ifdef REG_PARSE */
/*--------------------------------------------------------------------*/
#if 0

int reg_show (REGRESS *reg, FILE *file, int prec, double eps)
{                               /* --- describe a regression function */
  int    i, k, n;               /* loop variables, counter */
  double *c;                    /* to traverse the coefficients */
  int    *e;                    /* table of exponents */
  char   *s;                    /* separators between terms */
  int    w;                     /* width of number output */

  assert(reg && file            /* check the function arguments */
      && (prec >= 0) && (eps >= 0));
  if (prec <= 0) prec = 6;      /* adapt the number of places */
  if (reg->cnt *prec <= 40) {   /* if the output fits into one line, */
    s = " ";      w = 0; }      /* print only a blank between terms, */
  else {                        /* if it does not fit into one line, */
    s = "\n    "; w = prec+7; } /* print one term per line */
  if (reg->max > 0)             /* if logistic regression */
    fprintf(file, "z = %g /(1 +exp(y))\n", reg->max);
  fputs("y = ", file); n = 0;   /* start the function output */
  c = reg->cfs;                 /* traverse the reg. coefficients */
  if      (reg->dim <= 2) {     /* --- if single variable regression */
    for (i = reg->deg+1; --i > 1; c++) {
      if (fabs(*c) < eps) continue;   /* skip negligible terms */
      if (n++ > 0) fputs(s, file);    /* print a separator */
      fprintf(file, "%+*.*g x^%d", w, prec, *c, i);
    }                           /* print coeffs. of different powers */
    if (fabs(*c) >= eps) {      /* if linear term is not negligible */
      if (n++ > 0) fputs(s, file);
      fprintf(file, "%+*.*g x", w, prec, *c);
    } c++; }                    /* print the linear term */
  else if (reg->deg <= 1) {     /* --- if linear regression */
    for (i = 0; ++i < reg->dim; c++) {
      if (fabs(*c) < eps) continue;   /* skip negligible terms */
      if (n++ > 0) fputs(s, file);    /* print a separator */
      fprintf(file, "%+*.*g x_%d", w, prec, *c, i);
    } }                         /* print the linear terms */
  else if (reg->deg <= 2) {     /* --- if quadratic regression */
    for (i = 0; ++i < reg->dim; c++) {
      for (k = 0; ++k < i; c++) {
        if (fabs(*c) < eps) continue; /* skip negligible terms */
        if (n++ > 0) fputs(s, file);  /* print a separator */
        fprintf(file, "%+*.*g x_%d x_%d", w, prec, *c, k, i);
      }                         /* print the mixed terms */
      if (fabs(*c) < eps) continue;   /* skip negligible terms */
      if (n++ > 0) fputs(s, file);    /* print a separator */
      fprintf(file, "%+*.*g x_%d^2", w, prec, *c, i);
    }                           /* print the squares */
    for (i = 0; ++i < reg->dim; c++) {
      if (fabs(*c) < eps) continue;   /* skip negligible terms */
      if (n++ > 0) fputs(s, file);    /* print a separator */
      fprintf(file, "%+*.*g x_%d", w, prec, *c, i);
    } }                         /* print the linear terms */
  else {                        /* --- if general regression */
    e = reg->tab +(i = reg->dim-1);
    *--e = reg->deg +1;         /* init. the table of exponents to */
    while (--i > 0) *--e = 0;   /* x_n^(deg+1), i.e. "previous" state */
    for (i = reg->cnt; --i > 0; c++) {
      if (e[0] > 0) {           /* update the table of exponents: */
        e[1]++; e[0]--; }       /* if first exponent is not zero, */
      else {                    /* shift a 1 from first to second */
        for (k = 0; e[++k] <= 0; );   /* find first nonzero expo. */
        e[0] = e[k]-1; e[k] = 0; e[k+1]++;
      }                         /* move exponents back to first */
      if (fabs(*c) < eps) continue;   /* skip negligible terms */
      if (n++ > 0) fputs(s, file);    /* print a separator */
      fprintf(file, "%+*.*g", w, prec, *c);
      for (k = 0; ++k < reg->dim; ) {
        if      (e[k-1] == 1) fprintf(file, " x_%d", k);
        else if (e[k-1] >  1) fprintf(file, " x_%d^%d", k, e[k-1]);
      }                         /* print coefficients of */
    }                           /* the different products and */
  }                             /* a description the products */
  if (fabs(*c) >= eps) {        /* if constant term is not negligible */
    if (n++ > 0) fputs(s, file);/* print the constant term */
    fprintf(file, "%+*.*g", w, prec, *c); }
  else if (n <= 0)              /* if no term printed yet, */
    fputc('0', file);           /* print a zero as a fallback */
  fputc('\n', file);            /* terminate the output and */
  return ferror(file) ? -1 : 0; /* return the write status */
}  /* reg_show() */

#endif