classifier.cc

   1
   2 //////////////////////////////////////////////////////////////////////////////////
   3 // This program is free software: you can redistribute it and/or modify         //
   4 // it under the terms of the version 3 of the GNU General Public License        //
   5 // as published by the Free Software Foundation.                                //
   6 //                                                                              //
   7 // This program is distributed in the hope that it will be useful, but          //
   8 // WITHOUT ANY WARRANTY; without even the implied warranty of                   //
   9 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU             //
  10 // General Public License for more details.                                     //
  11 //                                                                              //
  12 // You should have received a copy of the GNU General Public License            //
  13 // along with this program. If not, see <http://www.gnu.org/licenses/>.         //
  14 //                                                                              //
  15 // Written by Francois Fleuret                                                  //
  16 // Copyright (C) Ecole Polytechnique Federale de Lausanne                       //
  17 // Contact <francois.fleuret@epfl.ch> for comments & bug reports                //
  18 //////////////////////////////////////////////////////////////////////////////////
  19
  20 // $Id: classifier.cc,v 1.3 2007-08-23 08:36:50 fleuret Exp $
  21
  22 #include <stdlib.h>
  23 #include <cmath>
  24
  25 #include "classifier.h"
  26
  27 DataSet::DataSet(const DataSet &ds) {
  28   copy(ds);
  29 }
  30
  31 DataSet::DataSet(int ns, int nf) {
  32   nb_samples = ns;
  33   nb_features = nf;
  34
  35   size = (nb_samples + 31)/32;
  36
  37   raw = new RawData;
  38
  39   raw->nrefs = 1;
  40   raw->x = new uint32_t[nb_features * size];
  41   raw->y = new uint32_t[size];
  42
  43   y_va = raw->y;
  44   x_va = new uint32_t *[nb_features];
  45   for(int j = 0; j < nb_features; j++) x_va[j] = raw->x + size*j;
  46 }
  47
  48 DataSet::DataSet(istream &is) {
  49   is >> nb_samples >> nb_features;
  50   size = (nb_samples + 31)/32;
  51
  52   raw = new RawData;
  53
  54   raw->nrefs = 1;
  55   raw->x = new uint32_t[nb_features * size];
  56   raw->y = new uint32_t[size];
  57
  58   y_va = raw->y;
  59   x_va = new uint32_t *[nb_features];
  60   for(int j = 0; j < nb_features; j++) x_va[j] = raw->x + size*j;
  61
  62   int v;
  63
  64   for(int s = 0; s < nb_samples; s++) {
  65     //    cerr << "s=" << s << "\n"; cerr.flush();
  66     for(int f = 0; f < nb_features; f++) {
  67       is >> v;
  68       fe_set_bit(s, x_va[f], v > 0);
  69     }
  70
  71     is >> v;
  72     fe_set_bit(s, y_va, v > 0);
  73
  74     if(is.eof()) {
  75       cerr << "Error: missing data!\n";
  76       exit(1);
  77     }
  78   }
  79 }
  80
  81 DataSet::DataSet(const DataSet &ds, const FeatureSelector &fs) {
  82   nb_samples = ds.nb_samples;
  83   nb_features = fs.nb_selected_features;
  84   size = (nb_samples + 31)/32;
  85
  86   raw = ds.raw;
  87
  88   raw->nrefs++;
  89
  90   y_va = raw->y;
  91   x_va = new uint32_t *[nb_features];
  92   for(int j = 0; j < nb_features; j++) x_va[j] = raw->x + size*fs.selected_index[j];
  93 }
  94
  95 DataSet::DataSet(const DataSet &ds, bool *selected_samples) {
  96   nb_samples = 0;
  97   for(int i = 0; i < ds.nb_samples; i++) if(selected_samples[i]) nb_samples++;
  98   nb_features = ds.nb_features;
  99   size = (nb_samples + 31)/32;
 100   raw = new RawData;
 101
 102   raw->nrefs = 1;
 103   raw->x = new uint32_t[nb_features * size];
 104   raw->y = new uint32_t[size];
 105
 106   y_va = raw->y;
 107   x_va = new uint32_t *[nb_features];
 108   for(int j = 0; j < nb_features; j++) x_va[j] = raw->x + size*j;
 109
 110   int k = 0;
 111   for(int s = 0; s < ds.nb_samples; s++) if(selected_samples[s]) {
 112     fe_set_bit(k, y_va, fe_get_bit(s, ds.y_va));
 113     for(int j = 0; j < nb_features; j++) fe_set_bit(k, x_va[j], fe_get_bit(s, ds.x_va[j]));
 114     k++;
 115   }
 116 }
 117
 118 DataSet &DataSet::operator = (const DataSet &ds) {
 119   if(this != &ds) {
 120     delete[] x_va;
 121     raw->nrefs--;
 122     if(raw->nrefs == 0) {
 123       delete[] raw->x;
 124       delete[] raw->y;
 125       delete raw;
 126     }
 127     copy(ds);
 128   }
 129   return *this;
 130 }
 131
 132 DataSet::~DataSet() {
 133   delete[] x_va;
 134   raw->nrefs--;
 135   if(raw->nrefs == 0) {
 136     delete[] raw->x;
 137     delete[] raw->y;
 138     delete raw;
 139   }
 140 }
 141
 142 void DataSet::copy(const DataSet &ds) {
 143   nb_samples = ds.nb_samples;
 144   nb_features = ds.nb_features;
 145   size = (nb_samples + 31)/32;
 146   raw = ds.raw;
 147   raw->nrefs++;
 148   y_va = raw->y;
 149   x_va = new uint32_t *[nb_features];
 150   for(int j = 0; j < nb_features; j++) x_va[j] = ds.x_va[j];
 151 }
 152
 153 void DataSet::save_ascii(ostream &os) {
 154   os << nb_samples << " " << nb_features << "\n";
 155   for(int s = 0; s < nb_samples; s++) {
 156     for(int f = 0; f < nb_features; f++) {
 157       if(fe_get_bit(s, x_va[f])) os << "1"; else os << "0";
 158       if(f < nb_features-1) os << " "; else os << "\n";
 159     }
 160     if(fe_get_bit(s, y_va)) os << "1\n"; else os << "0\n";
 161   }
 162 }
 163
 164 Classifier *Classifier::load(istream &is) {
 165   Classifier *result;
 166
 167   int scheme;
 168   is >> scheme;
 169
 170   switch(scheme) {
 171   case Classifier::ID_LINEAR:
 172     result = new LinearClassifier();
 173     break;
 174
 175   default:
 176     result = 0;
 177     cerr << "Unknown classifier type!\n";
 178     exit(1);
 179   }
 180
 181   result->inner_load(is);
 182
 183   return result;
 184 }
 185
 186 Classifier::~Classifier() {}
 187
 188 FeatureSelector::FeatureSelector(int nb) : nb_selected_features(nb),
 189                                            selected_index(new int[nb_selected_features]),
 190                                            weights(new float[nb_selected_features]) { }
 191
 192 FeatureSelector::~FeatureSelector() {
 193   delete[] weights;
 194   delete[] selected_index;
 195 }
 196
 197 LinearClassifier::LinearClassifier(int nf) : nb_features(nf),
 198                                              weights(new float[nb_features]),
 199                                              bias(0) { }
 200
 201 LinearClassifier::LinearClassifier() : nb_features(0),
 202                                        weights(0),
 203                                        bias(0) { }
 204
 205 LinearClassifier::~LinearClassifier() {
 206   delete[] weights;
 207 }
 208
 209 void LinearClassifier::compute_bayesian_weights(int nb_samples, uint32_t *y_va, uint32_t **x_va) {
 210   for(int nf = 0; nf < nb_features; nf++) {
 211     int n11 = fe_count_and(nb_samples, y_va, x_va[nf]);
 212     int n10 = fe_count_and_not(nb_samples, y_va, x_va[nf]);
 213     int n01 = fe_count_and_not(nb_samples, x_va[nf], y_va);
 214     int n00 = fe_count_and_not_not(nb_samples, y_va, x_va[nf]);
 215     if(n00 == 0) n00 = 1; // This is sort of a dirty way to emulate +/- infty
 216     if(n01 == 0) n01 = 1;
 217     if(n10 == 0) n10 = 1;
 218     if(n11 == 0) n11 = 1;
 219     weights[nf] = log(float(n11 * n00) / float(n10 * n01));
 220   }
 221 }
 222
 223 void LinearClassifier::compute_bias(int nb_samples, uint32_t *y_va, uint32_t **x_va, bool balanced_error) {
 224   Couple tmp[nb_samples];
 225
 226   int n00 = 0, n01 = 0, n10 = 0, n11 = 0;
 227   for(int s = 0; s < nb_samples; s++) {
 228     float r = 0;
 229     for(int nf = 0; nf < nb_features; nf++) if(fe_get_bit(s, x_va[nf])) r += weights[nf];
 230     tmp[s].index = s;
 231     tmp[s].value = r;
 232     if(fe_get_bit(s, y_va)) n11++; else n01++;
 233   }
 234
 235   qsort(tmp, nb_samples, sizeof(Couple), compare_couple);
 236
 237   float error, best_error = 2.0;
 238   bias = 0;
 239   for(int t = 0; t < nb_samples-1; t++) {
 240     if(fe_get_bit(tmp[t].index, y_va)) { n10++; n11--; } else { n01--; n00++; }
 241
 242     if(balanced_error)
 243       error = (float(n01) / float(n00 + n01) + float(n10) / float(n10 + n11)) / 2.0;
 244     else
 245       error = float(n01+n10) / float(n00 + n01 + n10 + n11);
 246
 247     if(error < best_error) {
 248       best_error = error;
 249       bias = - (tmp[t].value + tmp[t+1].value)/2.0;
 250     }
 251   }
 252 }
 253
 254 void LinearClassifier::learn_bayesian(const DataSet &ds, bool balanced_error) {
 255   compute_bayesian_weights(ds.nb_samples, ds.y_va, ds.x_va);
 256   compute_bias(ds.nb_samples, ds.y_va, ds.x_va, balanced_error);
 257 }
 258
 259 void LinearClassifier::learn_perceptron(const DataSet &ds, bool balanced_error) {
 260   for(int i = 0; i < nb_features; i++) weights[i] = 0.0;
 261
 262   int n_loop_max = 5000;
 263
 264   for(int i = 0; i < n_loop_max * ds.nb_samples; i++) {
 265     int ns = i % ds.nb_samples;
 266     float r = 0;
 267     for(int f = 0; f < nb_features; f++)
 268       if(fe_get_bit(ns, ds.x_va[f])) r += weights[f]; else r -= weights[f];
 269     float correct;
 270     if(fe_get_bit(ns, ds.y_va)) correct = 1.0; else correct = -1.0;
 271     if((r < 0 && correct >= 0) || (r >= 0 && correct < 0)) {
 272       for(int f = 0; f < nb_features; f++)
 273         if(fe_get_bit(ns, ds.x_va[f])) weights[f] += correct; else weights[f] += -correct;
 274     }
 275   }
 276
 277   compute_bias(ds.nb_samples, ds.y_va, ds.x_va, balanced_error);
 278 }
 279
 280 void FeatureSelector::cmim(const DataSet &ds) {
 281   fe_selection_cmim(ds.nb_samples, ds.nb_features, ds.x_va, ds.y_va, nb_selected_features, selected_index);
 282 }
 283
 284 void FeatureSelector::mim(const DataSet &ds) {
 285   fe_selection_mim(ds.nb_samples, ds.nb_features, ds.x_va, ds.y_va, nb_selected_features, selected_index);
 286 }
 287
 288 void FeatureSelector::random(const DataSet &ds) {
 289   bool used[ds.nb_features];
 290   for(int i = 0; i < ds.nb_features; i++) used[i] = false;
 291   int f;
 292   for(int nf = 0; nf < nb_selected_features; nf++) {
 293     do { f = int(drand48() * ds.nb_features); } while(used[f]);
 294     used[f] = true;
 295     selected_index[nf] = f;
 296   }
 297 }
 298
 299 FeatureSelector::FeatureSelector(istream &is) {
 300   is >> nb_selected_features;
 301   weights = new float[nb_selected_features];
 302   selected_index = new int[nb_selected_features];
 303   for(int i = 0; i < nb_selected_features; i++) is >> selected_index[i];
 304 }
 305
 306 void FeatureSelector::save(ostream &os) {
 307   os << nb_selected_features << "\n";
 308   for(int i = 0; i < nb_selected_features; i++) os << selected_index[i] << ((i < nb_selected_features-1) ? " " : "\n");
 309 }
 310
 311 void LinearClassifier::inner_load(istream &is) {
 312   is >> nb_features;
 313   delete[] weights;
 314   weights = new float[nb_features];
 315   for(int i = 0; i < nb_features; i++) is >> weights[i];
 316   is >> bias;
 317 }
 318
 319 void LinearClassifier::save(ostream &os) {
 320   os << ID_LINEAR << "\n";
 321   os << nb_features << "\n";
 322   for(int i = 0; i < nb_features; i++) os << weights[i] << ((i < nb_features-1) ? " " : "\n");
 323   os << bias << "\n";
 324 }
 325
 326 void LinearClassifier::predict(const DataSet &ds, float *result) {
 327   for(int s = 0; s < ds.nb_samples; s++) {
 328     float r = bias;
 329     for(int nf = 0; nf < nb_features; nf++) if(fe_get_bit(s, ds.x_va[nf])) r += weights[nf];
 330     result[s] = r;
 331   }
 332 }
 333
 334 void compute_error_rates(FeatureSelector *selector, Classifier *classifier,
 335                          const DataSet &testing_set, int &n00, int &n01, int &n10, int &n11, float *result) {
 336
 337   DataSet reduced_test_set(testing_set, *selector);
 338
 339   classifier->predict(reduced_test_set, result);
 340
 341   n00 = 0; n01 = 0; n10 = 0; n11 = 0;
 342   for(int s = 0; s < testing_set.nb_samples; s++) {
 343     if(fe_get_bit(s, testing_set.y_va)) {
 344       if(result[s] >= 0) n11++; else n10++;
 345     } else {
 346       if(result[s] >= 0) n01++; else n00++;
 347     }
 348   }
 349 }