6e2ecf0ca103a0b236b2d4af3786c13a0485425d
[cmim.git] / 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 }