v0.3.9b:
[mldemos:baraks-mldemos.git] / _AlgorithmsPlugins / ESMLR / classifierESMLR.cpp
1 /*********************************************************************
2 MLDemos: A User-Friendly visualization toolkit for machine learning
3 Copyright (C) 2010  Basilio Noris
4 Contact: mldemos@b4silio.com
5
6 Evolution-Strategy Mixture of Logisitics Regression
7 Copyright (C) 2011  Stephane Magnenat
8 Contact: stephane at magnenat dot net
9
10 This library is free software; you can redistribute it and/or
11 modify it under the terms of the GNU Lesser General Public License,
12 version 3 as published by the Free Software Foundation.
13
14 This library is distributed in the hope that it will be useful, but
15 WITHOUT ANY WARRANTY; without even the implied warranty of
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
17 Lesser General Public License for more details.
18
19 You should have received a copy of the GNU Lesser General Public
20 License along with this library; if not, write to the Free
21 Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
22 *********************************************************************/
23
24 #include "classifierESMLR.h"
25
26 #include <Eigen/Core>
27 #include <Eigen/Eigen>
28 #include <vector>
29 #include <map>
30 #include <iostream>
31
32 // ES-related stuff
33
34 namespace ESMLR
35 {
36         using namespace Eigen;
37         
38         double uniformRand(double min, double max)
39         {
40                 const double v = double(rand())/double(RAND_MAX);
41                 return (min + v * (max-min));
42         }
43
44         double gaussianRand(double mean, double sigm)
45         {
46                 // Generation using the Polar (Box-Mueller) method.
47                 // Code inspired by GSL, which is a really great math lib.
48                 // http://sources.redhat.com/gsl/
49                 // C++ wrapper available.
50                 // http://gslwrap.sourceforge.net/
51                 double r, x, y;
52
53                 // Generate random number in unity circle.
54                 do
55                 {
56                         x = uniformRand(-1, 1);
57                         y = uniformRand(-1, 1);
58                         r = x*x + y*y;
59                 }
60                 while (r > 1.0 || r == 0);
61
62                 // Box-Muller transform.
63                 return sigm * y * sqrt (-2.0 * log(r) / r) + mean;
64         }
65         
66         double sigm(double v)
67         {
68                 return 1. / (1. + exp(-v));
69         }
70         
71         double sgn(double v)
72         {
73                 if (v > 0) return 1;
74                 if (v < 0) return -1;
75                 return 0;
76         }
77         
78         struct Classifier
79         {
80                 MatrixXd w;
81                 VectorXd w_b;
82                 VectorXd v;
83                 double v_b;
84                 double beta;
85                 
86                 Classifier(unsigned cutCount, unsigned dataSize, double beta):
87                         w(cutCount, dataSize),
88                         w_b(cutCount),
89                         v(cutCount),
90                         beta(beta)
91                 {}
92                 double eval(const VectorXd& x) const;
93                 friend std::ostream& operator<< (std::ostream& stream, const Classifier& that);
94         };
95         
96         double Classifier::eval(const VectorXd& x) const
97         {
98                 assert(w.rows() == w_b.size());
99                 assert(w_b.size() == v.size());
100                 assert(w.cols() == x.size());
101                 double sum(0);
102                 for (int i = 0; i < w.rows(); ++i)
103                         sum += v(i) * (2. * sigm(beta * (w.row(i).dot(x) + w_b(i))) - 1.);
104                 sum += v_b;
105                 return 2. * sigm(w.rows() * 2 * sum) - 1.;
106                 //return sum > 0 ? 1 : -1;
107         }
108         
109         std::ostream& operator<< (std::ostream& stream, const Classifier& that)
110         {
111                 stream << "Classifier on " << that.w.cols() << " dimensions with " << that.w.rows() << " hyperplanes" << std::endl;
112                 stream << "w:\n" << that.w << std::endl;
113                 stream << "w_b:\n" << that.w_b << std::endl;
114                 stream << "v:\n" << that.v << std::endl;
115                 stream << "v_b:\n" << that.v_b << std::endl;
116                 stream << "beta:\n" << that.beta << std::endl;
117                 return stream;
118         }
119         
120         struct Individual
121         {
122                 Classifier classifier;
123                 
124                 double r_w;
125                 double r_w_b;
126                 double r_v;
127                 double r_v_b;
128                 
129                 Individual(unsigned cutCount = 0, unsigned dataSize = 0, double beta = 1);
130                 void mutate(double dataAvrSd);
131                 Individual createChild(double dataAvrSd) const;
132                 static Individual createRandom(unsigned cutCount, unsigned dataSize, double dataAvrSd, double beta);
133         };
134         
135         Individual::Individual(unsigned cutCount, unsigned dataSize, double beta):
136                 classifier(cutCount, dataSize, beta),
137                 r_w(1),
138                 r_w_b(1),
139                 r_v(1),
140                 r_v_b(1)
141         {}
142         
143         void Individual::mutate(double dataAvrSd)
144         {
145                 assert(classifier.w.rows() == classifier.w_b.size());
146                 assert(classifier.w.rows() == classifier.v.size());
147                 // mutate rate
148                 r_w *= (rand()%2 == 0) ? 1.25 : 0.8;
149                 r_w_b *= (rand()%2 == 0) ? 1.25 : 0.8;
150                 r_v *= (rand()%2 == 0) ? 1.25 : 0.8;
151                 r_v_b *= (rand()%2 == 0) ? 1.25 : 0.8;
152                 // mutate using rate
153                 const double aprioriRate(0.05);
154                 // w
155                 for (int i = 0; i < classifier.w.rows(); ++i)
156                 {
157                         for (int j = 0; j < classifier.w.cols(); ++j)
158                                 classifier.w(i,j) += gaussianRand(0, aprioriRate*r_w);
159                         classifier.w.row(i) /= classifier.w.row(i).norm();
160                 }
161                 // w_b
162                 for (int i = 0; i < classifier.w_b.size(); ++i)
163                         classifier.w_b(i) += gaussianRand(0, dataAvrSd*aprioriRate*r_w_b);
164                 // v
165                 for (int i = 0; i < classifier.v.size(); ++i)
166                         classifier.v(i) += gaussianRand(0, aprioriRate*r_v);
167                 classifier.v /= classifier.v.norm();
168                 // v_b
169                 classifier.v_b += gaussianRand(0, double(classifier.v.size())*aprioriRate*r_v_b);
170         }
171         
172         Individual Individual::createChild(double dataAvrSd) const
173         {
174                 Individual child(*this);
175                 child.mutate(dataAvrSd);
176                 return child;
177         }
178         
179         Individual Individual::createRandom(unsigned cutCount, unsigned dataSize, double dataAvrSd, double beta)
180         {
181                 Individual ind(cutCount, dataSize, beta);
182                 Classifier& classifier(ind.classifier);
183                 // w
184                 for (int i = 0; i < classifier.w.rows(); ++i)
185                 {
186                         for (int j = 0; j < classifier.w.cols(); ++j)
187                                 classifier.w(i,j) = uniformRand(-1, 1);
188                         classifier.w.row(i) /= classifier.w.row(i).norm();
189                 }
190                 // w_b
191                 for (int i = 0; i < classifier.w_b.size(); ++i)
192                         classifier.w_b(i) = gaussianRand(0, dataAvrSd);
193                 // v
194                 for (int i = 0; i < classifier.v.size(); ++i)
195                         classifier.v(i) = uniformRand(-1, 1);
196                 classifier.v /= classifier.v.norm();
197                 // v_b
198                 classifier.v_b = uniformRand(-double(classifier.v.size()), double(classifier.v.size()));
199                 return ind;
200         }
201         
202         struct Population: protected std::vector<Individual>
203         {
204                 typedef std::pair<double, double> ErrorPair;
205                 
206                 Population(unsigned cutCount, unsigned dataSize, double dataAvrSd, double beta, unsigned indPerDim);
207                 ErrorPair evolveOneGen(const VectorXd& y, const MatrixXd& x, double dataAvrSd);
208                 Classifier optimise(const VectorXd& y, const MatrixXd& x, double dataAvrSd, size_t genCount);
209         };
210         
211         Population::Population(unsigned cutCount, unsigned dataSize, double dataAvrSd, double beta, unsigned indPerDim):
212         std::vector<Individual>((cutCount*(dataSize+1)+1)*indPerDim)
213         {
214                 // create initial population
215                 for (iterator it(begin()); it != end(); ++it)
216                         *it = Individual::createRandom(cutCount, dataSize, dataAvrSd, beta);
217         }
218         
219         Population::ErrorPair Population::evolveOneGen(const VectorXd& y, const MatrixXd& x, double dataAvrSd)
220         {
221                 assert(y.size() == x.rows());
222                 typedef std::multimap<double, Individual> EvaluationMap;
223                 typedef EvaluationMap::iterator EvaluationMapIterator;
224                 EvaluationMap evalutationMap;
225                 
226                 // evaluation
227                 double totalError(0);
228                 for (const_iterator it(begin()); it != end(); ++it)
229                 {
230                         const Individual& ind(*it);
231                         double error(0);
232                         for (int sample = 0; sample < y.size(); ++sample)
233                         {
234                                 const double v(ind.classifier.eval(x.row(sample)));
235                                 const double delta(y(sample) - v);
236                                 error += delta*delta;
237                         }
238                         totalError += error;
239                         evalutationMap.insert(std::make_pair(error, ind));
240                 }
241                 const double averageError(totalError / double(size()));
242                 
243         // selection
244                 assert((size() / 4) * 4 == size());
245                 size_t ind = 0;
246                 for (EvaluationMapIterator it = evalutationMap.begin(); ind < size() / 4; ++it, ++ind)
247                 {
248                         //cout << "S " << it->first << "\n";
249                         (*this)[ind * 4] = it->second;
250                         (*this)[ind * 4 + 1] = it->second.createChild(dataAvrSd);
251                         (*this)[ind * 4 + 2] = it->second.createChild(dataAvrSd);
252                         (*this)[ind * 4 + 3] = it->second.createChild(dataAvrSd);
253                 }
254                 
255                 // return statistics
256                 return ErrorPair(evalutationMap.begin()->first, averageError);
257         }
258         
259         Classifier Population::optimise(const VectorXd& y, const MatrixXd& x, double dataAvrSd, size_t genCount)
260         {
261                 // optimise
262                 for (size_t g = 0; g < genCount; ++g)
263                 {
264                         const ErrorPair e = evolveOneGen(y, x, dataAvrSd);
265                         std::cout << g << " : " << e.first << ", " << e.second << ", ";
266                         // compute number of missclassified
267                         unsigned missClassified(0);
268                         for (int sample = 0; sample < y.size(); ++sample)
269                                 missClassified += fabs(sgn((*this)[0].classifier.eval(x.row(sample))) - y(sample)) / 2;
270                         std::cout << missClassified << std::endl;
271                 }
272                 return (*this)[0].classifier;
273         }
274 }
275
276 // classifier for mld
277
278 ClassifierESMLR::ClassifierESMLR():
279         cutCount(2),
280         alpha(3),
281         genCount(30),
282         indPerDim(50),
283         classifier(0)
284 {
285         bSingleClass = true;
286         //bUseDrawTimer = false;
287         classThresh = 0;
288         classSpan = 1;
289         dim = 0;
290 }
291
292 ClassifierESMLR::~ClassifierESMLR()
293 {
294         if (classifier)
295                 delete classifier;
296 }
297
298 void ClassifierESMLR::Train(std::vector< fvec > samples, ivec labels)
299 {
300         assert(labels.size() == samples.size());
301         dim = samples[0].size();
302         Eigen::VectorXd y(labels.size());
303         Eigen::MatrixXd x(samples.size(), dim);
304         for (size_t i = 0; i < labels.size(); ++i)
305         {
306                 y(i) = labels[i];
307                 for (size_t j = 0; j < samples[i].size(); ++j)
308                         x(i,j) = samples[i][j];
309         }
310         std::cerr << "Learning data:\ny:\n" << y << "\nx:\n" << x << std::endl;
311         
312         // compute stddev
313         const Eigen::VectorXd avr = x.colwise().sum() / double(x.rows());
314         Eigen::VectorXd stdDev(Eigen::VectorXd::Constant(avr.size(), 0));
315         for (int i = 0; i < x.rows(); ++i)
316         {
317                 const Eigen::VectorXd& row(x.row(i));
318                 const Eigen::VectorXd delta(row - avr);
319                 stdDev += delta.cwiseProduct(delta);
320         }
321         stdDev /= x.rows();
322         stdDev = stdDev.cwiseSqrt();
323         double dataAvrSd(stdDev.sum() / double(stdDev.size()));
324         std::cerr << "stddev:\n" << stdDev << "\navr stddev: " << dataAvrSd << std::endl;
325         
326         const double beta(alpha / dataAvrSd);
327         ESMLR::Population pop(cutCount, dim, dataAvrSd, beta, indPerDim);
328         if (classifier)
329                 delete classifier;
330         classifier = new ESMLR::Classifier(pop.optimise(y, x, dataAvrSd, genCount));
331         std::cerr << *classifier << std::endl;
332 }
333
334 float ClassifierESMLR::Test(const fvec &sample)
335 {
336         assert(classifier);
337         Eigen::VectorXd x(dim);
338         FOR(d, dim) x(d) = sample[d];
339         return classifier->eval(x);
340 }
341
342 char *ClassifierESMLR::GetInfoString()
343 {
344         char *text = new char[1024];
345         snprintf(text, 1024, "Evolution Strategy, Mixture of Logistic Regressions\n"
346         "hyperplane count: %d\n"
347         "alpha: %f\n"
348         "generation count: %d\n"
349         "individual per dim: %d\n",
350         cutCount, alpha, genCount, indPerDim
351         );
352         return text;
353 }
354
355 void ClassifierESMLR::SetParams(u32 cutCount, float alpha, u32 genCount, u32 indPerDim)
356 {
357         this->cutCount = cutCount;
358         this->alpha = alpha;
359         this->genCount = genCount;
360         this->indPerDim = indPerDim;
361 }