fastentropy.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: fastentropy.cc,v 1.4 2007-08-23 08:36:50 fleuret Exp $
  21
  22 using namespace std;
  23
  24 #include <iostream>
  25 #include <cmath>
  26 #include <stdio.h>
  27 #include <stdlib.h>
  28
  29 #include "misc.h"
  30 #include "fastentropy.h"
  31
  32 int fe_nb_bits[65536];
  33 double fe_logn[65536], fe_nlogn[65536];
  34
  35 #ifdef DEBUG
  36 int fe_was_initialized = 0;
  37 #endif
  38
  39 inline void fe_and(int n, uint32_t *a, uint32_t *b, uint32_t *ab) {
  40   for(int k = 0; k < (n+31)/32; k++) *ab++ = *a++ & *b++;
  41 }
  42
  43 inline void fe_and_not(int n, uint32_t *a, uint32_t *b, uint32_t *ab) {
  44   for(int k = 0; k < (n+31)/32; k++) *ab++ = *a++ & ~*b++;
  45 }
  46
  47 inline void fe_and_not_not(int n, uint32_t *a, uint32_t *b, uint32_t *ab) {
  48   for(int k = 0; k < (n+31)/32; k++) *ab++ = ~*a++ & ~*b++;
  49 }
  50
  51 inline int fe_count(int n, uint32_t *a) {
  52   uint16_t *aa = (uint16_t *) a;
  53   int t = 0;
  54   for(int k = 0; k < n/16; k++) t += fe_nb_bits[*aa++];
  55   if(n%16 > 0) t += fe_nb_bits[(65535 >> (16-(n%16))) & *aa];
  56   return t;
  57 }
  58
  59 inline double fe_entropy(int n, uint32_t *a) {
  60   int n1 = fe_count(n, a);
  61   return fe_logn[n] - (fe_nlogn[n1] + fe_nlogn[n-n1])/double(n);
  62 }
  63
  64 inline double fe_entropy_couple(int n, uint32_t *a, uint32_t *b) {
  65   int n11 = fe_count_and(n, a, b);
  66   int n10 = fe_count_and_not(n, a, b);
  67   int n01 = fe_count_and_not(n, b, a);
  68   int n00 = fe_count_and_not_not(n, a, b);
  69
  70   return fe_logn[n] - (fe_nlogn[n00] + fe_nlogn[n01] + fe_nlogn[n10] + fe_nlogn[n11])/double(n);
  71 }
  72
  73 void fe_init_tables() {
  74 #ifdef DEBUG
  75   fe_was_initialized = 1;
  76 #endif
  77   for(int i = 0; i < 65536; i++) {
  78     if(i == 0) { fe_logn[i] = 0.0; fe_nlogn[i] = 0.0; }
  79     else { fe_logn[i] = log(double(i)); fe_nlogn[i] = double(i) * log(double(i)); }
  80     int n = 0;
  81     for(int j = 0; j < 16; j++) if(i & (1 << j)) n++;
  82     fe_nb_bits[i] = n;
  83   }
  84 }
  85
  86 void fe_selection_cmim(int nb_samples,
  87                        int nb_tests, uint32_t **x, uint32_t *y,
  88                        int nb_selected, int *selected) {
  89
  90 #ifdef DEBUG
  91   if(!fe_was_initialized) {
  92     cerr << "fe_init_tables() was not called!\n";
  93     abort();
  94   }
  95 #endif
  96
  97   if(nb_samples > 65535) {
  98     cerr << "Too many pictures, the nlogn table is too small.\n";
  99     exit(1);
 100   }
 101
 102   double s[nb_tests];
 103   double ch[nb_tests];
 104   int m[nb_tests];
 105
 106   double h = fe_entropy(nb_samples, y);
 107
 108   for(int i = 0; i < nb_tests; i++) {
 109     ch[i] = fe_entropy_couple(nb_samples, y, x[i]) - fe_entropy(nb_samples, x[i]);
 110     s[i] = h - ch[i];
 111     m[i] = 0;
 112   }
 113
 114   for(int n = 0; n < nb_selected; n++) {
 115     double best_s = 0;
 116     for(int i = 0; i < nb_tests; i++) {
 117       if(s[i] > best_s) {
 118         int nb_uint32 = (nb_samples+31)/32;
 119         if(m[i] < n) {
 120           uint32_t z00[nb_uint32], z01[nb_uint32], z10[nb_uint32], z11[nb_uint32];
 121           fe_and(nb_samples, x[i], y, z11);
 122           fe_and_not(nb_samples, x[i], y, z10);
 123           fe_and_not(nb_samples, y, x[i], z01);
 124           fe_and_not_not(nb_samples, x[i], y, z00);
 125           while(s[i] > best_s && m[i] < n) {
 126             double h_y_xi_xmi = fe_logn[nb_samples] -
 127               (  fe_nlogn[fe_count_and    (nb_samples, z11, x[selected[m[i]]])]
 128                + fe_nlogn[fe_count_and    (nb_samples, z10, x[selected[m[i]]])]
 129                + fe_nlogn[fe_count_and    (nb_samples, z01, x[selected[m[i]]])]
 130                + fe_nlogn[fe_count_and    (nb_samples, z00, x[selected[m[i]]])]
 131                + fe_nlogn[fe_count_and_not(nb_samples, z11, x[selected[m[i]]])]
 132                + fe_nlogn[fe_count_and_not(nb_samples, z10, x[selected[m[i]]])]
 133                + fe_nlogn[fe_count_and_not(nb_samples, z01, x[selected[m[i]]])]
 134                + fe_nlogn[fe_count_and_not(nb_samples, z00, x[selected[m[i]]])])/double(nb_samples);
 135             double h_xi_xmi = fe_entropy_couple(nb_samples, x[i], x[selected[m[i]]]);
 136             double ss = ch[selected[m[i]]] - (h_y_xi_xmi - h_xi_xmi);
 137             if(ss < s[i]) s[i] = ss;
 138             m[i]++;
 139           }
 140         }
 141         if(s[i] > best_s) {
 142           best_s = s[i];
 143           selected[n] = i;
 144         }
 145       }
 146     }
 147   }
 148 }
 149
 150 void fe_selection_mim(int nb_samples,
 151                       int nb_tests, uint32_t **x, uint32_t *y,
 152                       int nb_selected, int *selected) {
 153
 154 #ifdef DEBUG
 155   if(!fe_was_initialized) {
 156     cerr << "fe_init_tables() was not called!\n";
 157     abort();
 158   }
 159 #endif
 160
 161   if(nb_samples > 65535) {
 162     cerr << "Too many pictures, the nlogn table is too small.\n";
 163     exit(1);
 164   }
 165
 166   Couple tmp[nb_tests];
 167   double h = fe_entropy(nb_samples, y);
 168
 169   for(int i = 0; i < nb_tests; i++) {
 170     tmp[i].index = i;
 171     tmp[i].value = fe_entropy_couple(nb_samples, y, x[i]) - h - fe_entropy(nb_samples, x[i]);
 172   }
 173
 174   qsort(tmp, nb_tests, sizeof(Couple), compare_couple);
 175
 176   // Here we have the features sorted according to their mutual
 177   // information with the class to predict
 178
 179   for(int n = 0; n < nb_selected; n++) selected[n] = tmp[n].index;
 180 }