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. //
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. //
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/>. //
15 // Written by Francois Fleuret //
16 // Copyright (C) Ecole Polytechnique Federale de Lausanne //
17 // Contact <francois.fleuret@epfl.ch> for comments & bug reports //
18 //////////////////////////////////////////////////////////////////////////////////
20 // $Id: fastentropy.cc,v 1.4 2007-08-23 08:36:50 fleuret Exp $
30 #include "fastentropy.h"
32 int fe_nb_bits[65536];
33 double fe_logn[65536], fe_nlogn[65536];
36 int fe_was_initialized = 0;
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++;
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++;
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++;
51 inline int fe_count(int n, uint32_t *a) {
52 uint16_t *aa = (uint16_t *) a;
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];
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);
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);
70 return fe_logn[n] - (fe_nlogn[n00] + fe_nlogn[n01] + fe_nlogn[n10] + fe_nlogn[n11])/double(n);
73 void fe_init_tables() {
75 fe_was_initialized = 1;
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)); }
81 for(int j = 0; j < 16; j++) if(i & (1 << j)) n++;
86 void fe_selection_cmim(int nb_samples,
87 int nb_tests, uint32_t **x, uint32_t *y,
88 int nb_selected, int *selected) {
91 if(!fe_was_initialized) {
92 cerr << "fe_init_tables() was not called!\n";
97 if(nb_samples > 65535) {
98 cerr << "Too many pictures, the nlogn table is too small.\n";
106 double h = fe_entropy(nb_samples, y);
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]);
114 for(int n = 0; n < nb_selected; n++) {
116 for(int i = 0; i < nb_tests; i++) {
118 int nb_uint32 = (nb_samples+31)/32;
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;
150 void fe_selection_mim(int nb_samples,
151 int nb_tests, uint32_t **x, uint32_t *y,
152 int nb_selected, int *selected) {
155 if(!fe_was_initialized) {
156 cerr << "fe_init_tables() was not called!\n";
161 if(nb_samples > 65535) {
162 cerr << "Too many pictures, the nlogn table is too small.\n";
166 Couple tmp[nb_tests];
167 double h = fe_entropy(nb_samples, y);
169 for(int i = 0; i < nb_tests; i++) {
171 tmp[i].value = fe_entropy_couple(nb_samples, y, x[i]) - h - fe_entropy(nb_samples, x[i]);
174 qsort(tmp, nb_tests, sizeof(Couple), compare_couple);
176 // Here we have the features sorted according to their mutual
177 // information with the class to predict
179 for(int n = 0; n < nb_selected; n++) selected[n] = tmp[n].index;