mpr.c

it is regression Algorithm in C/C++.

        switch (*s++) {         /* evaluate option */
          case 'M': matinp = 1;                     break;
          case 'o': optarg = &trgname;              break;
          case 'x': deg    = (int)strtol(s, &s, 0); break;
          case 'l': logit  =      strtod(s, &s);    break;
          case 'e': mode  |= REG_EXPOS;             break;
          case 'v': verb   = 1;                     break;
          case 'b': optarg = &blanks;               break;
          case 'f': optarg = &fldseps;              break;
          case 'r': optarg = &recseps;              break;
          case 'C': optarg = &comment;              break;
          case 'd': flags |= AS_DFLT;               break;
          case 'h': optarg = &fn_hdr;               break;
          default : error(E_OPTION, *--s);          break;
        }                       /* set option variables */
        if (!*s) break;         /* if at end of string, abort loop */
        if (optarg) { *optarg = s; optarg = NULL; break; }
      } }                       /* get option argument */
    else {                      /* -- if argument is no option */
      switch (k++) {            /* evaluate non-option */
        case  0: fn_dom = s;      break;
        case  1: fn_tab = s;      break;
        case  2: fn_reg = s;      break;
        default: error(E_ARGCNT); break;
      }                         /* note file names */
    }
  }
  if (optarg) error(E_OPTARG);  /* check option argument */
  if (matinp) {                 /* if matrix version */
    if (k != 2) error(E_ARGCNT);/* check the number of arguments */
    fn_reg = fn_tab;            /* shift the file names */
    fn_tab = fn_dom; fn_dom = NULL; }
  else {                        /* if table version */
    if (k != 3) error(E_ARGCNT);/* check number of arguments */
    if (fn_hdr && (strcmp(fn_hdr, "-") == 0))
      fn_hdr = "";              /* convert "-" to "" */
    i = (!fn_dom || !*fn_dom) ? 1 : 0;
    if  (!fn_tab || !*fn_tab) i++;
    if  ( fn_hdr && !*fn_hdr) i++;
    if (i > 1) error(E_STDIN);  /* stdin must not be used twice */
    if (fn_hdr) flags = AS_ATT | (flags & ~AS_DFLT);
  }                             /* set the header file flag */
  if (deg <= 0) error(E_DEGREE, deg);
  
  if (matinp) {                 /* if matrix version */
    /* --- create a table file scanner --- */
    t     = clock();            /* start the timer */
    tscan = ts_create();        /* create a table file scanner */
    if (!tscan) error(E_NOMEM); /* set the separator characters */
    if (blanks)  ts_chars(tscan, TS_BLANK,   blanks);
    if (fldseps) ts_chars(tscan, TS_FLDSEP,  fldseps);
    if (recseps) ts_chars(tscan, TS_RECSEP,  recseps);
    if (comment) ts_chars(tscan, TS_COMMENT, comment);
    ts_chars(tscan,TS_NULL,""); /* remove the null value characters */
    err = ts_info(tscan);       /* get the error information */

    /* --- read the matrix file --- */
    if (fn_tab && *fn_tab)      /* if a file name is given, */
      in = fopen(fn_tab, "r");  /* open the file for reading */
    else {                      /* if no file name is given, */
      in = stdin; fn_tab = "<stdin>"; }   /* read from stdin */
    fprintf(stderr, "\nreading %s ... ", fn_tab);
    if (!in) error(E_FOPEN, fn_tab);
    vec = vec_readx(tscan, in, &dim);
    if (!vec) {                 /* read the first vector */
      if (err->code >= 0) error(E_FREAD, fn_tab);
      error(err->code, fn_tab, 1, err->s, err->fld, err->exp);
    }
    if (dim < 2) error(E_DIM);  /* check the number of dimensions */
    reg = reg_create(dim, deg); /* (there must be at least 2) */
    if (!reg) error(E_NOMEM);   /* create a regression object */
    if (logit > 0)              /* if logistic regression, */
      reg_setmax(reg, logit);   /* set the transformation parameter */
    reg_init(reg);              /* initialize the regression object */
    do {                        /* aggregate and count a vector */
      if (reg_aggr(reg, vec, 1) != 0)
        error(E_LOGIT, vec[dim-1]);
      if (((++cnt & 0x03ff) == 0) && verb)
        fprintf(stderr, "%8d\b\b\b\b\b\b\b\b", cnt);
    } while (vec_read(vec, dim, tscan, in) == 0);
    if (err->code < 0)          /* check for a read error */
      error(err->code, fn_tab, cnt+1, err->s, err->fld, err->exp);
    fprintf(stderr, "[%d vector(s)] ", cnt);
    fprintf(stderr, "done [%.2fs].\n", SEC_SINCE(t)); }

  else {                        /* if table version */
    /* --- read attribute set --- */
    t    = clock();             /* start the timer */
    scan = sc_create(fn_dom);   /* create a scanner */
    if (!scan) error((!fn_dom || !*fn_dom) ? E_NOMEM : E_FOPEN, fn_dom);
    attset = as_create("domains", att_delete);
    if (!attset) error(E_NOMEM);/* create an attribute set */
    fprintf(stderr, "\nreading %s ... ", sc_fname(scan));
    if ((sc_nexter(scan)   <  0)/* start scanning (get first token) */
    ||  (as_parse(attset, scan, AT_ALL) != 0)
    ||  (as_attcnt(attset) <= 0)/* parse attribute set and */
    ||  !sc_eof(scan))          /* check for end of file */
      error(E_PARSE, sc_fname(scan));
    attcnt = as_attcnt(attset); /* print a success message */
    fprintf(stderr, "[%d attribute(s)] ", attcnt);
    fprintf(stderr, "done [%.2fs].", SEC_SINCE(t));

    /* --- filter attributes --- */
    for (i = attcnt; --i >= 0; ) {
      att = as_att(attset, i);  /* traverse the attributes */
      k   = att_getdir(att);    /* and get their directions */
      att_setmark(att, ((k == DIR_IN) || (k == DIR_OUT)) ? -1 : +1);
    }                           /* mark input and output attributes */
    as_attcut(NULL, attset, AS_MARKED);
    attcnt = as_attcnt(attset); /* cut all other attributes and */
    for (i = attcnt; --i >= 0;) /* clear all attribute markers */
      att_setmark(as_att(attset, i), 0);

    /* --- determine id of target column --- */
    if (trgname) {              /* if a target att. name is given */
      for (i = attcnt; --i >= 0; ) /* remove attribute directions */
        att_setdir(as_att(attset, i), DIR_IN);
      trgid = as_attid(attset, trgname);
      if (trgid < 0) error(E_TARGET, trgname, sc_fname(scan)); }
    else {                      /* if no target att. name is given */
      for (trgid = -1, i = attcnt; --i >= 0; ) {
        if (att_getdir(as_att(attset, i)) != DIR_OUT) continue;
        if (trgid < 0) trgid = i; else error(E_MULTTRG);
      }                         /* find a (unique) output attribute */
      if (trgid < 0) trgid = attcnt -1;
    }                           /* by default use the last attribute */
    att = as_att(attset,trgid); /* get the target attribute and */
    att_setdir(att, DIR_OUT);   /* set the target attribute direction */
    sc_delete(scan); scan = NULL;  /* delete the scanner */

    /* --- create regression object --- */
    if (as_attcnt(attset) < 2)  /* check the number of variables */
      error(E_DIM);             /* (there must be at least 2) */
    attmap = am_create(attset, 0, 1.0);
    if (!attmap) error(E_NOMEM);/* create an attribute map */
    am_target(attmap, trgid);   /* set and check */
    if (am_outcnt(attmap) != 1) /* the target attribute */
      error(E_TARGET, att_name(as_att(attset, trgid)));
    reg = reg_createx(attmap, deg);
    if (!reg) error(E_NOMEM);   /* create a regression object */
    if (logit > 0)              /* if logistic regression, */
      reg_setmax(reg, logit);   /* set the transformation parameter */
    reg_init(reg);              /* initialize the regression object */
    fprintf(stderr, "\n");      /* terminate the startup message */

    /* --- read table header --- */
    t  = clock();               /* start the timer */
    as_chars(attset, recseps, fldseps, blanks, "", comment);
    in = io_hdr(attset, fn_hdr, fn_tab, flags|AS_MARKED, 1);
    if (!in) error(1);          /* read the table header */

    /* --- process tuples --- */
    f = AS_INST | (flags & ~(AS_ATT|AS_DFLT));
    i = ((flags & AS_DFLT) && !(flags & AS_ATT))
      ? 0 : as_read(attset, in, f);
    while (i == 0) {            /* record read loop */
      tplcnt++;                 /* count tuple and sum its weight */
      tplwgt += as_getwgt(attset);
      if (reg_aggrx(reg, NULL) != 0)
        error(E_LOGIT, (double)att_inst(am_att(attmap, -1))->f);
      if (((++cnt & 0x03ff) == 0) && verb)
        fprintf(stderr, "%8d\b\b\b\b\b\b\b\b", cnt);
      i = as_read(attset,in,f); /* aggregate the data tuple and */
    }                           /* try to read the next record */
    if (i < 0) {                /* if an error occurred, */
      err = as_err(attset);     /* get the error information */
      tplcnt += (flags & (AS_ATT|AS_DFLT)) ? 1 : 2;
      io_error(i, fn_tab, tplcnt, err->s, err->fld, err->exp);
      error(1);                 /* print an error message */
    }                           /* and abort the program */
    fprintf(stderr, "[%d/%g tuple(s)] ", tplcnt, tplwgt);
    fprintf(stderr, "done [%.2fs].\n", SEC_SINCE(t));
  }                             /* print a success message */

  /* --- compute regression function --- */
  t = clock();                  /* start the timer */
  fprintf(stderr, "computing regression ... ");
  if (reg_solve(reg) < 0)       /* compute regression coefficients */
    error(E_REGRESS);           /* by solving the equation system */
  fprintf(stderr, "[sse: %g] ", reg_sse(reg));
  fprintf(stderr, "done [%.2fs].\n", SEC_SINCE(t));

  /* --- write the regression model --- */
  t = clock();                  /* start the timer */
  if (fn_reg && *fn_reg)        /* if an output file name is given, */
    out = fopen(fn_reg, "w");   /* open the output file */
  else {                        /* if no output file name is given, */
    out = stdout; fn_reg = "<stdout>"; }    /* write to std. output */
  fprintf(stderr, "writing %s ... ", fn_reg);
  if (!out) error(E_FOPEN, fn_reg);
  if (!matinp) {                /* if table version */
    if (as_desc(attset, out, AS_TITLE|AS_MARKED|AS_IVALS, 0) != 0)
      error(E_FWRITE, fn_reg);  /* describe attribute domains */
    fprintf(out, "\n");         /* leave one line empty */
  }
  if (reg_desc(reg, out, mode, 0) != 0)
    error(E_FWRITE, fn_reg);    /* describe the regression polynomial */
  if (out != stdout) {          /* if not written to standard output, */
    i = fclose(out); out = NULL;/* close the output file */
    if (i != 0) error(E_FWRITE, fn_reg);
  }                             /* print a success message */
  fprintf(stderr, "[%d parameters] ", reg_cnt(reg));
  fprintf(stderr, "done [%.2fs].\n", SEC_SINCE(t));

  /* --- clean up --- */
  #ifndef NDEBUG
  if (!matinp)                  /* if table version */
    reg_deletex(reg, 1);        /* delete the regression object */
  else {                        /* if matrix version */
    reg_delete(reg);            /* delete the regression object, */
    free(vec);                  /* the data vector buffer, */
    ts_delete(tscan);           /* and the table scanner */
  }                             /* if (matinp) .. else .. */
  #endif
  #ifdef STORAGE
  showmem("at end of program"); /* check memory usage */
  #endif
  return 0;                     /* return `ok' */
}  /* main() */