data-tool.cc

   1
   2 /*
   3  *  data-tool is a command line tool to do simple statistical
   4  *  processing on numerical data.
   5  *
   6  *  Copyright (c) 2009 Francois Fleuret
   7  *  Written by Francois Fleuret <francois@fleuret.org>
   8  *
   9  *  This file is part of data-tool.
  10  *
  11  *  data-tool is free software: you can redistribute it and/or modify
  12  *  it under the terms of the GNU General Public License version 3 as
  13  *  published by the Free Software Foundation.
  14  *
  15  *  data-tool is distributed in the hope that it will be useful, but
  16  *  WITHOUT ANY WARRANTY; without even the implied warranty of
  17  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  18  *  General Public License for more details.
  19  *
  20  *  You should have received a copy of the GNU General Public License
  21  *  along with data-tool.  If not, see <http://www.gnu.org/licenses/>.
  22  *
  23  */
  24
  25 #include <iostream>
  26 #include <cmath>
  27 #include <stdlib.h>
  28 #include <string.h>
  29
  30 using namespace std;
  31
  32 struct Couple {
  33   int index;
  34   double value;
  35 };
  36
  37 int compare_couple(const void *a, const void *b) {
  38   if(((Couple *) a)->value < ((Couple *) b)->value) return -1;
  39   else if(((Couple *) a)->value > ((Couple *) b)->value) return 1;
  40   else return 0;
  41 }
  42
  43 double *inflate_array(double *x, int current_size, int new_size) {
  44   double *xx = new double[new_size];
  45   for(int n = 0; n < current_size; n++) xx[n] = x[n];
  46   delete[] x;
  47   return xx;
  48 }
  49
  50 char *next_word(char *buffer, char *r, int buffer_size) {
  51   char *s;
  52   s = buffer;
  53   if(r != NULL)
  54     {
  55       if(*r == '"') {
  56         r++;
  57         while((*r != '"') && (*r != '\0') &&
  58               (s<buffer+buffer_size-1))
  59           *s++ = *r++;
  60         if(*r == '"') r++;
  61       } else {
  62         while((*r != '\r') && (*r != '\n') && (*r != '\0') &&
  63               (*r != '\t') && (*r != ' ') && (*r != ',') &&
  64               (s<buffer+buffer_size-1))
  65           *s++ = *r++;
  66       }
  67
  68       while((*r == ' ') || (*r == '\t') || (*r == ',')) r++;
  69       if((*r == '\0') || (*r=='\r') || (*r=='\n')) r = NULL;
  70     }
  71   *s = '\0';
  72   return r;
  73 }
  74
  75 void check_opt(int argc, char **argv, int n_opt, int n, const char *help) {
  76   if(n_opt + n >= argc) {
  77     cerr << "ERROR: Missing argument for " << argv[n_opt] << ". Expecting " << help << "." << endl;
  78     exit(1);
  79   }
  80 }
  81
  82 void print_help_and_exit(int e) {
  83   cout << "Simple data processing tool. Written by Francois Fleuret." << endl
  84        << endl
  85        << "This application takes data from the standard input and prints" << endl
  86        << "the result on the standard output. It expects either a list of" << endl
  87        << "float values (to produce histograms, cumulative distribution functions" << endl
  88        << "or the mean, variance, etc.) or a list of couples of values of the form" << endl
  89        << "x y on each line (where the sign of x tells the class and y the parameter" << endl
  90        << "value) to compute the ROC curve or the ROC curve surface.\n" << endl
  91        << "The options are:" << endl
  92        << "  --help" << endl
  93        << "  --roc" << endl
  94        << "  --roc-surface" << endl
  95        << "  --error" << endl
  96        << "  --normalize" << endl
  97        << "  --histo" << endl
  98        << "  --cumul" << endl
  99        << "  --misc" << endl
 100        << "  --auto-extrema" << endl
 101        << "  --xbounds <float: xmin> <float: xmax>" << endl
 102        << "  --ybounds <float: ymin> <float: ymax>" << endl
 103        << "  --nb-bins <int: number of bins>" << endl;
 104   exit(e);
 105 }
 106
 107 void check_single_processing(bool unknown_processing) {
 108   if(!unknown_processing) {
 109     cerr << "ERROR: You can't do two different processings." << endl;
 110     exit(1);
 111   }
 112 }
 113
 114 int main(int argc, char **argv) {
 115   double xmin = 0, xmax = 1, ymin = 0, ymax = 1;
 116   int nb_bins = 10;
 117   const int buffer_size = 1024;
 118
 119   char line[buffer_size], token[buffer_size];
 120   bool auto_extrema = false;
 121   bool normalize = false;
 122
 123   int i = 1;
 124
 125   enum { UNKNOWN, ROC, ROC_SURFACE, ERROR, HISTO, CUMUL, MISC } processing = UNKNOWN;
 126
 127   // Parsing the command line arguments ////////////////////////////////
 128
 129   while(i < argc) {
 130
 131     if(argc == 1 || strcmp(argv[i], "--help") == 0) print_help_and_exit(0);
 132
 133     else if(strcmp(argv[i], "--roc") == 0) {
 134       check_single_processing(processing == UNKNOWN);
 135       processing = ROC;
 136       i++;
 137     }
 138
 139     else if(strcmp(argv[i], "--roc-surface") == 0) {
 140       check_single_processing(processing == UNKNOWN);
 141       processing = ROC_SURFACE;
 142       i++;
 143     }
 144
 145     else if(strcmp(argv[i], "--error") == 0) {
 146       check_single_processing(processing == UNKNOWN);
 147       processing = ERROR;
 148       i++;
 149     }
 150
 151     else if(strcmp(argv[i], "--cumul") == 0) {
 152       check_single_processing(processing == UNKNOWN);
 153       processing = CUMUL;
 154       i++;
 155     }
 156
 157     else if(strcmp(argv[i], "--normalize") == 0) {
 158       normalize = true;
 159       i++;
 160     }
 161
 162     else if(strcmp(argv[i], "--histo") == 0) {
 163       check_single_processing(processing == UNKNOWN);
 164       processing = HISTO;
 165       i++;
 166     }
 167
 168     else if(strcmp(argv[i], "--misc") == 0) {
 169       check_single_processing(processing == UNKNOWN);
 170       processing = MISC;
 171       i++;
 172     }
 173
 174     else if(strcmp(argv[i], "--auto-extrema") == 0) {
 175       auto_extrema = true;
 176       i++;
 177     }
 178
 179     else if(strcmp(argv[i], "--xbounds") == 0) {
 180       check_opt(argc, argv, i, 2, "<float: xmin> <float: xmax>");
 181       xmin = atof(argv[i+1]);
 182       xmax = atof(argv[i+2]);
 183       if(xmin >= xmax) {
 184         cerr << "ERROR: Incorrect bounds." << endl;
 185         exit(1);
 186       }
 187       i += 3;
 188     }
 189
 190     else if(strcmp(argv[i], "--ybounds") == 0) {
 191       check_opt(argc, argv, i, 2, "<float: ymin> <float: ymax>");
 192       ymin = atof(argv[i+1]);
 193       ymax = atof(argv[i+2]);
 194       if(ymin >= ymax) {
 195         cerr << "ERROR: Incorrect bounds." << endl;
 196         exit(1);
 197       }
 198       i += 3;
 199     }
 200
 201     else if(strcmp(argv[i], "--nb-bins") == 0) {
 202       check_opt(argc, argv, i, 1, "<int: number of bins>");
 203       nb_bins = atoi(argv[i+1]);
 204       if(nb_bins < 1) {
 205         cerr << "ERROR: Incorrect number of bins." << endl;
 206         exit(1);
 207       }
 208       i += 2;
 209     }
 210
 211     else {
 212       cerr << "ERROR: Unknown option " << argv[i]  << endl;
 213       print_help_and_exit(1);
 214     }
 215   }
 216
 217   // Processing the data ///////////////////////////////////////////////
 218
 219   switch(processing) {
 220
 221   case CUMUL:
 222
 223     {
 224       int nb_samples = 0, nb_samples_max = 50000;
 225       double *x = new double[nb_samples_max];
 226
 227       while(!cin.eof()) {
 228         if(nb_samples == nb_samples_max) {
 229           x = inflate_array(x, nb_samples_max, 2 * nb_samples_max);
 230           nb_samples_max = 2 * nb_samples_max;
 231         }
 232
 233         cin.getline(line, buffer_size);
 234
 235         if(line[0]) {
 236           char *s = line;
 237           s = next_word(token, s, buffer_size);
 238           x[nb_samples] = atof(token);
 239           nb_samples++;
 240         }
 241       }
 242
 243       Couple *tmp = new Couple[nb_samples];
 244       for(int n = 0; n < nb_samples; n++) {
 245         tmp[n].index = n;
 246         tmp[n].value = x[n];
 247       }
 248
 249       qsort(tmp, nb_samples, sizeof(Couple), compare_couple);
 250
 251       for(int n = 0; n < nb_samples; n++)
 252         cout << tmp[n].value << " " << double(n+1)/double(nb_samples)  << endl;
 253
 254       delete[] tmp;
 255       delete[] x;
 256
 257     }
 258
 259     break;
 260
 261   case ROC:
 262   case ROC_SURFACE:
 263   case ERROR:
 264
 265     {
 266       int nb_samples = 0, nb_samples_max = 1000;
 267       double *x = new double[nb_samples_max], *y = new double[nb_samples_max];
 268
 269       while(!cin.eof()) {
 270         if(nb_samples == nb_samples_max) {
 271           x = inflate_array(x, nb_samples_max, 2 * nb_samples_max);
 272           y = inflate_array(y, nb_samples_max, 2 * nb_samples_max);
 273           nb_samples_max = 2 * nb_samples_max;
 274         }
 275
 276         cin.getline(line, buffer_size);
 277
 278         if(line[0]) {
 279           char *s = line;
 280           s = next_word(token, s, buffer_size);
 281           x[nb_samples] = atof(token);
 282           s = next_word(token, s, buffer_size);
 283           y[nb_samples] = atof(token);
 284           nb_samples++;
 285         }
 286       }
 287
 288       Couple *tmp = new Couple[nb_samples];
 289       int nb_rn = 0, nb_rp = 0, nb_fp = 0, nb_fn = 0;
 290
 291       bool binary = true;
 292       for(int n = 0; binary && n < nb_samples; n++) binary &= (x[n] == 0 || x[n] == 1);
 293       if(binary) {
 294         cerr << "WARNING: your classes are binary, I process them accordingly." << endl;
 295         for(int n = 0; n < nb_samples; n++) x[n] = 2 * x[n] - 1;
 296       }
 297
 298       for(int n = 0; n < nb_samples; n++) {
 299         tmp[n].index = n;
 300         tmp[n].value = y[n];
 301         if(x[n] >= 0) nb_rp++;
 302         else { nb_rn++; nb_fp++; }
 303       }
 304
 305       if(nb_rp == 0) cerr << "WARNING: No true positive." << endl;
 306       if(nb_rn == 0) cerr << "WARNING: No true negative." << endl;
 307
 308       qsort(tmp, nb_samples, sizeof(Couple), compare_couple);
 309
 310       if(processing == ROC) {
 311         for(int n = 0; n < nb_samples - 1; n++) {
 312           if(x[tmp[n].index] >= 0) nb_fn++;
 313           else                     nb_fp--;
 314           if(tmp[n].value < tmp[n+1].value) {
 315             cout << double(nb_fp)/double(nb_rn) << " "
 316                  << 1 - double(nb_fn) / double(nb_rp) << " "
 317                  << (tmp[n].value + tmp[n+1].value)/2 << " "
 318                  << endl;
 319           }
 320         }
 321       } else if(processing == ROC_SURFACE) {
 322         double surface = 0;
 323         double cx = double(nb_fp)/double(nb_rn), cy = 1 - double(nb_fn) / double(nb_rp);
 324         for(int n = 0; n < nb_samples - 1; n++) {
 325           if(x[tmp[n].index] >= 0) nb_fn++;
 326           else                     nb_fp--;
 327           if(tmp[n].value < tmp[n+1].value) {
 328             double ncx = double(nb_fp)/double(nb_rn), ncy = 1 - double(nb_fn) / double(nb_rp);
 329             surface += (cx - ncx) * cy;
 330             cx = ncx; cy = ncy;
 331           }
 332         }
 333         cout << surface  << endl;
 334       } else {
 335         for(int n = 0; n < nb_samples - 1; n++) {
 336           if(x[tmp[n].index] >= 0) nb_fn++;
 337           else                     nb_fp--;
 338           if(tmp[n].value < tmp[n+1].value) {
 339             cout << (tmp[n].value + tmp[n+1].value)/2 << " "
 340                  << double(nb_fp + nb_fn)/double(nb_rn + nb_rp) << " "
 341                  << endl;
 342           }
 343         }
 344       }
 345
 346       delete[] tmp;
 347       delete[] x;
 348       delete[] y;
 349
 350     }
 351
 352     break;
 353
 354   case HISTO:
 355
 356     {
 357       int nb_samples = 0, nb_samples_max = 1000;
 358       double *x = new double[nb_samples_max];
 359
 360       while(!cin.eof()) {
 361         if(nb_samples == nb_samples_max) {
 362           x = inflate_array(x, nb_samples_max, 2 * nb_samples_max);
 363           nb_samples_max = 2 * nb_samples_max;
 364         }
 365
 366         cin.getline(line, buffer_size);
 367
 368         if(line[0]) {
 369           char *s = line;
 370           s = next_word(token, s, buffer_size);
 371           x[nb_samples] = atof(token);
 372           if(auto_extrema) {
 373             if(nb_samples == 0 || x[nb_samples] > xmax) xmax = x[nb_samples];
 374             if(nb_samples == 0 || x[nb_samples] < xmin) xmin = x[nb_samples];
 375           }
 376           nb_samples++;
 377         }
 378       }
 379
 380       int *nb = new int[nb_bins];
 381       for(int n = 0; n < nb_bins; n++) nb[n] = 0;
 382
 383       int nb_total = 0;
 384       for(int s = 0; s < nb_samples; s++) {
 385         int n = int((x[s] - xmin)/(xmax - xmin) * nb_bins);
 386         if(n >= 0 && n < nb_bins) nb[n]++;
 387         else {
 388           cerr << "WARNING: value " << x[s] << " is out of histogram." << endl;
 389         }
 390         nb_total++;
 391       }
 392
 393       if(normalize) {
 394         for(int n = 0; n < nb_bins; n++)
 395           cout << xmin + ((xmax - xmin) * n) / double(nb_bins) << " "
 396                << (nb[n] / double(nb_total))/((xmax - xmin) / double(nb_bins))  << endl;
 397       } else {
 398         for(int n = 0; n < nb_bins; n++)
 399           cout << xmin + ((xmax - xmin) * n) / double(nb_bins) << " "
 400                << nb[n] / double(nb_total)  << endl;
 401       }
 402
 403       delete[] nb;
 404     }
 405
 406     break;
 407
 408   case MISC:
 409
 410     {
 411       int nb_samples = 0, nb_samples_max = 1000;
 412       double *x = new double[nb_samples_max];
 413       int nb = 0;
 414       double min = 0, max = 0;
 415       double sum = 0, sumsq = 0;
 416
 417       while(!cin.eof()) {
 418         if(nb_samples == nb_samples_max) {
 419           x = inflate_array(x, nb_samples_max, 2 * nb_samples_max);
 420           nb_samples_max = 2 * nb_samples_max;
 421         }
 422
 423         cin.getline(line, buffer_size);
 424         char *s = line;
 425         if(line[0]) {
 426           s = next_word(token, s, buffer_size);
 427           x[nb_samples] = atof(token);
 428           nb_samples++;
 429           double x = atof(token);
 430           if(nb == 0 || x > max) max = x;
 431           if(nb == 0 || x < min) min = x;
 432           sum += x;
 433           sumsq += x*x;
 434           nb++;
 435         }
 436       }
 437
 438       Couple *tmp = new Couple[nb_samples];
 439       for(int n = 0; n < nb_samples; n++) {
 440         tmp[n].index = n;
 441         tmp[n].value = x[n];
 442       }
 443
 444       qsort(tmp, nb_samples, sizeof(Couple), compare_couple);
 445
 446       delete[] x;
 447
 448       double mu = sum / double(nb);
 449       double sigma = (sumsq - sum * mu) / double(nb - 1);
 450       double stdd = sqrt(sigma);
 451
 452       cout << "MIN " << min
 453            << " MAX " << max
 454            << " MU " << mu
 455            << " SIGMA " << sigma
 456            << " STDD " << stdd
 457            << " SUM " << sum
 458            << " MEDIAN " << tmp[nb_samples/2].value
 459            << " QUANTILE0.1 " << tmp[int(nb_samples * 0.1)].value
 460            << " QUANTILE0.9 " << tmp[int(nb_samples * 0.9)].value
 461            << endl;
 462
 463       delete[] tmp;
 464     }
 465
 466     break;
 467
 468   default:
 469     cerr << "ERROR: You must choose a processing type." << endl;
 470     exit(1);
 471   }
 472
 473 }