const-fixing functions
[mldemos:mldemos.git] / _AlgorithmsPlugins / Projections / classifierLinear.cpp
1 /*********************************************************************\r
2 MLDemos: A User-Friendly visualization toolkit for machine learning\r
3 Copyright (C) 2010  Basilio Noris\r
4 Contact: mldemos@b4silio.com\r
5 \r
6 This library is free software; you can redistribute it and/or\r
7 modify it under the terms of the GNU Lesser General Public\r
8 License as published by the Free Software Foundation; either\r
9 version 2.1 of the License, or (at your option) any later version.\r
10 \r
11 This library is distributed in the hope that it will be useful,\r
12 but WITHOUT ANY WARRANTY; without even the implied warranty of\r
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU\r
14 Library General Public License for more details.\r
15 \r
16 You should have received a copy of the GNU Lesser General Public\r
17 License along with this library; if not, write to the Free\r
18 Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.\r
19 *********************************************************************/\r
20 #include "public.h"\r
21 #include "basicMath.h"\r
22 #include "classifierLinear.h"\r
23 #include "JnS/Matutil.h"\r
24 #include "JnS/JnS.h"\r
25 #include <roc.h>\r
26 #include <QDebug>\r
27 \r
28 using namespace std;\r
29 \r
30 void ClassifierLinear::Train( std::vector< fvec > samples, ivec labels )\r
31 {\r
32         if(!samples.size()) return;\r
33 \r
34         switch(linearType)\r
35         {\r
36         case 0:\r
37                 TrainPCA(samples, labels);\r
38                 break;\r
39     case 1:\r
40         TrainLDA(samples, labels, -1);\r
41         break;\r
42     case 2:\r
43         TrainLDA(samples, labels, 1);\r
44         break;\r
45     case 3:\r
46         TrainLDA(samples, labels, 0);\r
47                 break;\r
48         }\r
49 \r
50         vector<fvec> projected; projected.resize(samples.size());\r
51         FOR(i, samples.size())\r
52         {\r
53                 projected[i] = Project(samples[i]);\r
54         }\r
55 \r
56         int dim = samples[0].size();\r
57         meanPos.resize(dim,0);\r
58         meanNeg.resize(dim,0);\r
59         int cntPos=0, cntNeg=0;\r
60         FOR(i, samples.size())\r
61         {\r
62                 if(labels[i]==1)\r
63                 {\r
64                         FOR(d,dim) meanPos[d] += projected[i][d];\r
65                         cntPos++;\r
66                 }\r
67                 else\r
68                 {\r
69                         FOR(d,dim) meanNeg[d] += projected[i][d];\r
70                         cntNeg++;\r
71                 }\r
72         }\r
73         FOR(d,dim)\r
74         {\r
75                 if(cntPos) meanPos[d] /= cntPos;\r
76                 if(cntNeg) meanNeg[d] /= cntNeg;\r
77         }\r
78     bUsesDrawTimer = true;\r
79         minResponse = FLT_MAX;\r
80         maxResponse = -FLT_MAX;\r
81     midResponse = 0.5;\r
82         float minResp = FLT_MAX;\r
83         float maxResp = -FLT_MAX;\r
84         FOR(i, samples.size())\r
85         {\r
86                 float response = Test(samples[i]);\r
87                 if(minResp > response) minResp = response;\r
88                 if(maxResp < response) maxResp = response;\r
89         }\r
90         if(minResp == maxResp)\r
91         {\r
92                         maxResp = minResp + 0.5f;\r
93                         minResp -= 0.5f;\r
94         }\r
95         minResponse = minResp;\r
96         maxResponse = maxResp;\r
97 \r
98     if(linearType < 3)\r
99     {\r
100         vector<f32pair> roc(samples.size());\r
101         FOR(i, samples.size())\r
102         {\r
103             fVec point(samples[i][0] - meanAll[0], samples[i][1] - meanAll[1]);\r
104             float estimate = W*point;\r
105             float response = -(estimate - threshold);\r
106             if(minResponse != FLT_MAX)\r
107             {\r
108                 response = (response-minResponse)/fabs(maxResponse-minResponse); // 0-1 range\r
109             }\r
110             roc[i] = make_pair(response, labels[i]);\r
111         }\r
112         midResponse = GetBestThreshold(roc);\r
113     }\r
114 }\r
115 \r
116 float ClassifierLinear::Test(const fvec &sample ) const\r
117 {\r
118         float response = 0;\r
119     if(linearType < 4) // pca, lda, fisher\r
120         {\r
121         fVec point(sample[0] - meanAll.at(0), sample[1] - meanAll.at(1));\r
122                 float estimate = W*point;\r
123                 response = -(estimate - threshold);\r
124         }\r
125         else\r
126         {\r
127                 float distPos = 0, distNeg = 0;\r
128                 if(meanPos.size() == sample.size() && meanNeg.size() == sample.size())\r
129                 {\r
130                         fvec projected = Project(sample);\r
131 \r
132                         FOR(d,sample.size())\r
133                         {\r
134                 distPos += fabs(projected[d] - meanPos.at(d));\r
135                 distNeg += fabs(projected[d] - meanNeg.at(d));\r
136                         }\r
137                 }\r
138                 response = (distNeg - distPos);\r
139         }\r
140         if(minResponse != FLT_MAX)\r
141         {\r
142                 response = (response-minResponse)/fabs(maxResponse-minResponse); // 0-1 range\r
143         response = (response-midResponse)*6.f;\r
144         }\r
145     return response;\r
146 }\r
147 \r
148 const char *ClassifierLinear::GetInfoString() const\r
149 {\r
150         char *text = new char[1024];\r
151         sprintf(text, "");\r
152         switch(linearType)\r
153         {\r
154         case 0:\r
155                 sprintf(text, "%sPCA\n", text);\r
156                 break;\r
157     case 1:\r
158         sprintf(text, "%sMeansOnly\n", text);\r
159         break;\r
160     case 2:\r
161         sprintf(text, "%sLDA\n", text);\r
162         break;\r
163     case 3:\r
164                 sprintf(text, "%sFisher LDA\n", text);\r
165                 break;\r
166         default:\r
167                 sprintf(text, "%sNaive Bayes\n", text);\r
168                 break;\r
169         }\r
170     if(linearType < 4)\r
171         {\r
172                 sprintf(text, "%sProjection Direction:\n\t%.3f %.3f\n", text, W.x, W.y);\r
173         }\r
174         return text;\r
175 }\r
176 \r
177 void Invert(double *sigma, double *invSigma)\r
178 {\r
179         double det = sigma[0]*sigma[3] - sigma[1]*sigma[2];\r
180         if(det == 0) return; // non inversible\r
181         invSigma[0] = sigma[3] / det;\r
182         invSigma[1] = -sigma[1] / det;\r
183         invSigma[2] = -sigma[2] / det;\r
184         invSigma[3] = sigma[0] / det;\r
185 }\r
186 \r
187 fvec ClassifierLinear::Project(const fvec &sample) const\r
188 {\r
189         fvec newSample = sample;\r
190     if(linearType < 4) // pca, lda, fisher\r
191         {\r
192         fVec mean(meanAll.at(0), meanAll.at(1));\r
193                 fVec point(sample[0], sample[1]);\r
194                 float dot = W*(point-mean);\r
195                 fVec proj(dot*W.x, dot*W.y);\r
196                 //proj += cvVec2(.5f,.5f);\r
197                 proj += mean;\r
198                 newSample[0] = proj.x;\r
199                 newSample[1] = proj.y;\r
200         }\r
201         return newSample;\r
202 }\r
203 \r
204 void ClassifierLinear::SetParams(const u32 linearType )\r
205 {\r
206         this->linearType = linearType;\r
207         if(linearType == 1 || linearType == 2) bSingleClass = false;\r
208         else bSingleClass = true;\r
209 }\r
210 \r
211 void ClassifierLinear::TrainPCA(std::vector< fvec > samples, const ivec &labels)\r
212 {\r
213         u32 dim = 2;\r
214 \r
215         meanAll.resize(dim,0);\r
216         float **sigma = NULL;\r
217 \r
218         FOR(i, samples.size())\r
219         {\r
220                 meanAll += samples[i];\r
221         }\r
222         meanAll /= samples.size();\r
223 \r
224         fvec mean;\r
225         mean.resize(dim,0);\r
226 \r
227         // we want zero mean samples\r
228         FOR(i, samples.size()) samples[i] -= meanAll;\r
229 \r
230         GetCovariance(samples, mean, &sigma);\r
231 \r
232         float invSigma1[2][2] = {{sigma[1][1],-sigma[0][1]},{-sigma[0][1],sigma[0][0]}};\r
233         FOR(i,2)\r
234         {\r
235                 FOR(j,2)\r
236                 {\r
237                         invSigma1[i][j] /= sigma[0][0]*sigma[1][1] - sigma[0][1]*sigma[1][0];\r
238                 }\r
239         }\r
240 \r
241         float determinant = (invSigma1[1][1]+invSigma1[0][0])*(invSigma1[1][1]+invSigma1[0][0])-4.0f*(invSigma1[0][0]*invSigma1[1][1]-invSigma1[1][0]*invSigma1[0][1]);\r
242         float eigenvalue1, eigenvalue2;\r
243         if (determinant > 0){\r
244                 eigenvalue1 = (invSigma1[1][1] + invSigma1[0][0] + sqrtf(determinant)) / 2.0f;\r
245                 eigenvalue2 = (invSigma1[1][1] + invSigma1[0][0] - sqrtf(determinant)) / 2.0f;\r
246         }else{\r
247                 printf("determinant is not positive during calculation of eigenvalues !!");\r
248                 return;\r
249         }\r
250         fVec e1, e2, tmp;\r
251 \r
252         if(invSigma1[0][0] - eigenvalue1 != 0)\r
253                 e1.x = -invSigma1[0][1]/(invSigma1[0][0]-eigenvalue1);\r
254         else e1.x = 0;\r
255         e1.y = 1.0;\r
256 \r
257         if(invSigma1[0][0] - eigenvalue2 != 0)\r
258                 e2.x = -invSigma1[0][1]/(invSigma1[0][0]-eigenvalue2);\r
259         else e2.x = 0;\r
260         e2.y = 1.0;\r
261 \r
262         // swap the eigens\r
263         if (eigenvalue1 < eigenvalue2)\r
264         {\r
265                 fVec tmp = e1;\r
266                 e1 = e2;\r
267                 e2 = tmp;\r
268                 float lambda = eigenvalue1;\r
269                 eigenvalue1 = eigenvalue2;\r
270                 eigenvalue2 = lambda;\r
271         }\r
272 \r
273         e1=e1.normalize();\r
274         e2=e2.normalize();\r
275 \r
276         W = e2;\r
277         if(W.x < 0) W *= -1;\r
278 \r
279         W = W.normalize();\r
280 \r
281         KILL(sigma);\r
282 \r
283     threshold = 0;\r
284         bool inverted = false;\r
285         // we do the actual classification\r
286         u32 steps = 1000;\r
287         u32 minError = samples.size();\r
288         for( int c=0; c<steps; c++)\r
289         {\r
290                 u32 error = 0;\r
291                 u32 negError = 0;\r
292                 float thresh = 1.f / steps * c;\r
293                 FOR(i, samples.size())\r
294                 {\r
295                         fVec point(samples[i][0], samples[i][1]);\r
296                         float estimate = W*point;\r
297                         if(labels[i])\r
298                         {\r
299                                 if (estimate < thresh) error++;\r
300                                 else negError++;\r
301                         }\r
302                         else\r
303                         {\r
304                                 if (estimate >= thresh) error++;\r
305                                 else negError++;\r
306                         }\r
307                 }\r
308                 if(minError > min(negError, error))\r
309                 {\r
310                         inverted = negError > error;\r
311                         minError = min(negError, error);\r
312                         threshold = thresh;\r
313                 }\r
314         }\r
315 \r
316         float error = minError / (f32)samples.size();\r
317 }\r
318 \r
319 void ClassifierLinear::TrainLDA(std::vector< fvec > samples, const ivec &labels, int LDAType)\r
320 {\r
321         // we reduce the problem to a one vs many classification\r
322         u32 dim = 2;\r
323 \r
324         meanAll.resize(dim,0);\r
325         FOR(i, samples.size()) meanAll += samples[i];\r
326         meanAll /= samples.size();\r
327         FOR(i, samples.size()) samples[i] -= meanAll;\r
328 \r
329         vector<fvec> positives, negatives;\r
330         FOR(i, samples.size())\r
331         {\r
332                 if(labels[i]==1) positives.push_back(samples[i]);\r
333                 else negatives.push_back(samples[i]);\r
334         }\r
335 \r
336         fvec mean1, mean2;\r
337         mean1.resize(dim,0);\r
338         mean2.resize(dim,0);\r
339         float **sigma1 = NULL;\r
340         float **sigma2 = NULL;\r
341 \r
342         FOR(i, positives.size())\r
343         {\r
344                 FOR(j,dim)\r
345                 {\r
346                         mean1[j] += positives[i][j];\r
347                 }\r
348         }\r
349         mean1 /= (float)positives.size();\r
350 \r
351         FOR(i, negatives.size())\r
352         {\r
353                 FOR(j,dim)\r
354                 {\r
355                         mean2[j] += negatives[i][j];\r
356                 }\r
357         }\r
358         mean2 /= (float)negatives.size();\r
359 \r
360     if(LDAType==-1) // means only\r
361     {\r
362         W = fVec(mean2[0] - mean1[0],mean2[1] - mean1[1]);\r
363         W = W.normalize();\r
364         return;\r
365     }\r
366     else if(LDAType==0) // fisher LDA\r
367         {\r
368                 GetCovariance(positives, mean1, &sigma1);\r
369                 GetCovariance(negatives, mean2, &sigma2);\r
370         }\r
371     else // standard LDA\r
372         {\r
373                 vector<fvec> posnegs;\r
374                 FOR(i, positives.size()) posnegs.push_back(positives[i] - mean1);\r
375                 FOR(i, negatives.size()) posnegs.push_back(negatives[i] - mean2);\r
376                 fvec zero;\r
377                 zero.resize(dim,0);\r
378                 GetCovariance(posnegs, zero, &sigma1);\r
379                 sigma2 = sigma1;\r
380         }\r
381 \r
382         float sigma[2][2] = {{sigma1[0][0]+sigma2[0][0], sigma1[0][1]+sigma2[0][1]},{sigma1[1][0]+sigma2[1][0], sigma1[1][1]+sigma2[1][1]}};\r
383         float invSigma[2][2] = {{sigma[1][1],-sigma[1][0]},{-sigma[0][1],sigma[0][0]}};\r
384         FOR(i,2)\r
385         {\r
386                 FOR(j,2)\r
387                 {\r
388                         invSigma[i][j] /= sigma[0][0]*sigma[1][1] - sigma[0][1]*sigma[1][0];\r
389                 }\r
390         }\r
391 \r
392         float dM[2] = {mean2[0] - mean1[0],mean2[1] - mean1[1]};\r
393         float w[2] = { invSigma[0][0]*dM[0] + invSigma[0][1]*dM[1] , invSigma[1][0]*dM[0] + invSigma[1][1]*dM[1] };\r
394 \r
395         float n = sqrtf(w[0]*w[0] + w[1]*w[1]);\r
396         w[0] /= n; w[1] /= n;\r
397 \r
398 //      float c = w[0]*(mean1[0]+mean2[0])/2 + w[0]*(mean1[1]+mean2[1])/2;\r
399         W = fVec(w[0], w[1]);\r
400         W = W.normalize();\r
401 \r
402         if(sigma1 == sigma2)\r
403         {\r
404                 KILL(sigma1);\r
405         }\r
406         else\r
407         {\r
408                 KILL(sigma1);\r
409                 KILL(sigma2);\r
410         }\r
411 }\r
412 \r
413 void ClassifierLinear::GetCovariance(const vector<fvec> &samples, const fvec &mean, float ***covar)\r
414 {\r
415         int dim = mean.size();\r
416         float **cov = (*covar);\r
417         if(!cov)\r
418         {\r
419                 cov = new float*[dim];\r
420                 FOR(i, dim) cov[i] = new float[dim];\r
421                 (*covar) = cov;\r
422         }\r
423         FOR(i, dim)\r
424         {\r
425                 FOR(j,dim)\r
426                 {\r
427                         cov[i][j] = 0;\r
428                 }\r
429         }\r
430 \r
431         for (u32 i=0;i<samples.size();i++){\r
432                 float dX = samples[i][0] - mean[0];\r
433                 float dY = samples[i][1] - mean[1];\r
434                 cov[0][0] += dX*dX;\r
435                 cov[1][1] += dY*dY;\r
436                 cov[0][1] += dX*dY;\r
437         }\r
438         cov[0][0] /= (float)samples.size();\r
439         cov[1][1] /= (float)samples.size();\r
440         cov[0][1] /= (float)samples.size();\r
441         cov[1][0] = cov[0][1];\r
442 }\r
443 \r
444 ClassifierLinear::~ClassifierLinear()\r
445 {\r
446     if(Transf) free(Transf);\r
447 }\r
448 \r
449 void ClassifierLinear::TrainICA(std::vector< fvec > samples, const ivec &labels )\r
450 {\r
451     if(!samples.size()) return;\r
452     u32 dim = samples[0].size();\r
453         meanAll.resize(dim,0);\r
454         FOR(i, samples.size())\r
455         {\r
456                 meanAll += samples[i];\r
457         }\r
458         meanAll /= samples.size();\r
459 \r
460         const int nbsensors = dim;\r
461         const int nbsamples = samples.size(); \r
462 \r
463         if(!Transf)\r
464         {\r
465                 if ((Transf = (double *) calloc(nbsensors*nbsensors, sizeof(double))) == NULL) OutOfMemory() ;\r
466         }\r
467         double *Data, *Mixing, *Global;\r
468         if ((Data   = (double *) calloc(nbsensors*nbsamples, sizeof(double))) == NULL) OutOfMemory() ;\r
469         if ((Mixing = (double *) calloc(nbsensors*nbsensors, sizeof(double))) == NULL) OutOfMemory() ;\r
470         if ((Global = (double *) calloc(nbsensors*nbsensors, sizeof(double))) == NULL) OutOfMemory() ;\r
471 \r
472         FOR(i, samples.size())\r
473         {\r
474         FOR(d, nbsensors)\r
475         {\r
476             Data[i*nbsensors + d] = samples[i][d] - meanAll[d];\r
477             //          Data[i*nbsensors + d] = rand()/(float)RAND_MAX/5.;\r
478         }\r
479         }\r
480 \r
481     Identity(Mixing, nbsensors);\r
482     Mixing[0] = 2.0 ;\r
483         Jade(Transf, Data, nbsensors, nbsamples ) ; \r
484         //Shibbs(Transf, Data, nbsensors, nbsamples ) ;\r
485 \r
486         FOR(i,nbsensors*nbsensors) Transf[i] /= 10;\r
487 \r
488     projected = vector<fvec>(samples.size());\r
489     FOR(i, samples.size())\r
490     {\r
491         projected[i].resize(dim);\r
492         FOR(d, dim)\r
493         {\r
494             projected[i][d] = Data[i*nbsensors + d];\r
495         }\r
496     }\r
497 \r
498         free(Data);\r
499         //      free(Transf);  \r
500         free(Mixing);\r
501         free(Global);\r
502 \r
503     W = fVec(Transf[0], Transf[1*nbsensors]);\r
504 }\r