1 /*M///////////////////////////////////////////////////////////////////////////////////////
2 //
3 // IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
4 //
5 // By downloading, copying, installing or using the software you agree to this license.
6 // If you do not agree to this license, do not download, install,
7 // copy or use the software.
8 //
9 //
10 // License Agreement
11 // For Open Source Computer Vision Library
12 //
13 // Copyright (C) 2000, Intel Corporation, all rights reserved.
14 // Copyright (C) 2014, Itseez Inc, all rights reserved.
15 // Third party copyrights are property of their respective owners.
16 //
17 // Redistribution and use in source and binary forms, with or without modification,
18 // are permitted provided that the following conditions are met:
19 //
20 // * Redistribution's of source code must retain the above copyright notice,
21 // this list of conditions and the following disclaimer.
22 //
23 // * Redistribution's in binary form must reproduce the above copyright notice,
24 // this list of conditions and the following disclaimer in the documentation
25 // and/or other materials provided with the distribution.
26 //
27 // * The name of the copyright holders may not be used to endorse or promote products
28 // derived from this software without specific prior written permission.
29 //
30 // This software is provided by the copyright holders and contributors "as is" and
31 // any express or implied warranties, including, but not limited to, the implied
32 // warranties of merchantability and fitness for a particular purpose are disclaimed.
33 // In no event shall the Intel Corporation or contributors be liable for any direct,
34 // indirect, incidental, special, exemplary, or consequential damages
35 // (including, but not limited to, procurement of substitute goods or services;
36 // loss of use, data, or profits; or business interruption) however caused
37 // and on any theory of liability, whether in contract, strict liability,
38 // or tort (including negligence or otherwise) arising in any way out of
39 // the use of this software, even if advised of the possibility of such damage.
40 //
41 //M*/
42
43 #include "precomp.hpp"
44
45 #include <stdarg.h>
46 #include <ctype.h>
47
48 /****************************************************************************************\
49 COPYRIGHT NOTICE
50 ----------------
51
52 The code has been derived from libsvm library (version 2.6)
53 (http://www.csie.ntu.edu.tw/~cjlin/libsvm).
54
55 Here is the original copyright:
56 ------------------------------------------------------------------------------------------
57 Copyright (c) 2000-2003 Chih-Chung Chang and Chih-Jen Lin
58 All rights reserved.
59
60 Redistribution and use in source and binary forms, with or without
61 modification, are permitted provided that the following conditions
62 are met:
63
64 1. Redistributions of source code must retain the above copyright
65 notice, this list of conditions and the following disclaimer.
66
67 2. Redistributions in binary form must reproduce the above copyright
68 notice, this list of conditions and the following disclaimer in the
69 documentation and/or other materials provided with the distribution.
70
71 3. Neither name of copyright holders nor the names of its contributors
72 may be used to endorse or promote products derived from this software
73 without specific prior written permission.
74
75
76 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
77 ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
78 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
79 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR
80 CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
81 EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
82 PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
83 PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
84 LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
85 NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
86 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
87 \****************************************************************************************/
88
89 namespace cv { namespace ml {
90
91 typedef float Qfloat;
92 const int QFLOAT_TYPE = DataDepth<Qfloat>::value;
93
94 // Param Grid
checkParamGrid(const ParamGrid & pg)95 static void checkParamGrid(const ParamGrid& pg)
96 {
97 if( pg.minVal > pg.maxVal )
98 CV_Error( CV_StsBadArg, "Lower bound of the grid must be less then the upper one" );
99 if( pg.minVal < DBL_EPSILON )
100 CV_Error( CV_StsBadArg, "Lower bound of the grid must be positive" );
101 if( pg.logStep < 1. + FLT_EPSILON )
102 CV_Error( CV_StsBadArg, "Grid step must greater than 1" );
103 }
104
105 // SVM training parameters
106 struct SvmParams
107 {
108 int svmType;
109 int kernelType;
110 double gamma;
111 double coef0;
112 double degree;
113 double C;
114 double nu;
115 double p;
116 Mat classWeights;
117 TermCriteria termCrit;
118
SvmParamscv::ml::SvmParams119 SvmParams()
120 {
121 svmType = SVM::C_SVC;
122 kernelType = SVM::RBF;
123 degree = 0;
124 gamma = 1;
125 coef0 = 0;
126 C = 1;
127 nu = 0;
128 p = 0;
129 termCrit = TermCriteria( CV_TERMCRIT_ITER+CV_TERMCRIT_EPS, 1000, FLT_EPSILON );
130 }
131
SvmParamscv::ml::SvmParams132 SvmParams( int _svmType, int _kernelType,
133 double _degree, double _gamma, double _coef0,
134 double _Con, double _nu, double _p,
135 const Mat& _classWeights, TermCriteria _termCrit )
136 {
137 svmType = _svmType;
138 kernelType = _kernelType;
139 degree = _degree;
140 gamma = _gamma;
141 coef0 = _coef0;
142 C = _Con;
143 nu = _nu;
144 p = _p;
145 classWeights = _classWeights;
146 termCrit = _termCrit;
147 }
148
149 };
150
151 /////////////////////////////////////// SVM kernel ///////////////////////////////////////
152 class SVMKernelImpl CV_FINAL : public SVM::Kernel
153 {
154 public:
SVMKernelImpl(const SvmParams & _params=SvmParams ())155 SVMKernelImpl( const SvmParams& _params = SvmParams() )
156 {
157 params = _params;
158 }
159
getType() const160 int getType() const CV_OVERRIDE
161 {
162 return params.kernelType;
163 }
164
calc_non_rbf_base(int vcount,int var_count,const float * vecs,const float * another,Qfloat * results,double alpha,double beta)165 void calc_non_rbf_base( int vcount, int var_count, const float* vecs,
166 const float* another, Qfloat* results,
167 double alpha, double beta )
168 {
169 int j, k;
170 for( j = 0; j < vcount; j++ )
171 {
172 const float* sample = &vecs[j*var_count];
173 double s = 0;
174 for( k = 0; k <= var_count - 4; k += 4 )
175 s += sample[k]*another[k] + sample[k+1]*another[k+1] +
176 sample[k+2]*another[k+2] + sample[k+3]*another[k+3];
177 for( ; k < var_count; k++ )
178 s += sample[k]*another[k];
179 results[j] = (Qfloat)(s*alpha + beta);
180 }
181 }
182
calc_linear(int vcount,int var_count,const float * vecs,const float * another,Qfloat * results)183 void calc_linear( int vcount, int var_count, const float* vecs,
184 const float* another, Qfloat* results )
185 {
186 calc_non_rbf_base( vcount, var_count, vecs, another, results, 1, 0 );
187 }
188
calc_poly(int vcount,int var_count,const float * vecs,const float * another,Qfloat * results)189 void calc_poly( int vcount, int var_count, const float* vecs,
190 const float* another, Qfloat* results )
191 {
192 Mat R( 1, vcount, QFLOAT_TYPE, results );
193 calc_non_rbf_base( vcount, var_count, vecs, another, results, params.gamma, params.coef0 );
194 if( vcount > 0 )
195 pow( R, params.degree, R );
196 }
197
calc_sigmoid(int vcount,int var_count,const float * vecs,const float * another,Qfloat * results)198 void calc_sigmoid( int vcount, int var_count, const float* vecs,
199 const float* another, Qfloat* results )
200 {
201 int j;
202 calc_non_rbf_base( vcount, var_count, vecs, another, results,
203 2*params.gamma, 2*params.coef0 );
204 // TODO: speedup this
205 for( j = 0; j < vcount; j++ )
206 {
207 Qfloat t = results[j];
208 Qfloat e = std::exp(std::abs(t)); // Inf value is possible here
209 Qfloat r = (Qfloat)((e - 1.) / (e + 1.)); // NaN value is possible here (Inf/Inf or similar)
210 if (cvIsNaN(r))
211 r = std::numeric_limits<Qfloat>::infinity();
212 if (t < 0)
213 r = -r;
214 CV_DbgAssert(!cvIsNaN(r));
215 results[j] = r;
216 }
217 }
218
calc_rbf(int vcount,int var_count,const float * vecs,const float * another,Qfloat * results)219 void calc_rbf( int vcount, int var_count, const float* vecs,
220 const float* another, Qfloat* results )
221 {
222 double gamma = -params.gamma;
223 int j, k;
224
225 for( j = 0; j < vcount; j++ )
226 {
227 const float* sample = &vecs[j*var_count];
228 double s = 0;
229
230 for( k = 0; k <= var_count - 4; k += 4 )
231 {
232 double t0 = sample[k] - another[k];
233 double t1 = sample[k+1] - another[k+1];
234
235 s += t0*t0 + t1*t1;
236
237 t0 = sample[k+2] - another[k+2];
238 t1 = sample[k+3] - another[k+3];
239
240 s += t0*t0 + t1*t1;
241 }
242
243 for( ; k < var_count; k++ )
244 {
245 double t0 = sample[k] - another[k];
246 s += t0*t0;
247 }
248 results[j] = (Qfloat)(s*gamma);
249 }
250
251 if( vcount > 0 )
252 {
253 Mat R( 1, vcount, QFLOAT_TYPE, results );
254 exp( R, R );
255 }
256 }
257
258 /// Histogram intersection kernel
calc_intersec(int vcount,int var_count,const float * vecs,const float * another,Qfloat * results)259 void calc_intersec( int vcount, int var_count, const float* vecs,
260 const float* another, Qfloat* results )
261 {
262 int j, k;
263 for( j = 0; j < vcount; j++ )
264 {
265 const float* sample = &vecs[j*var_count];
266 double s = 0;
267 for( k = 0; k <= var_count - 4; k += 4 )
268 s += std::min(sample[k],another[k]) + std::min(sample[k+1],another[k+1]) +
269 std::min(sample[k+2],another[k+2]) + std::min(sample[k+3],another[k+3]);
270 for( ; k < var_count; k++ )
271 s += std::min(sample[k],another[k]);
272 results[j] = (Qfloat)(s);
273 }
274 }
275
276 /// Exponential chi2 kernel
calc_chi2(int vcount,int var_count,const float * vecs,const float * another,Qfloat * results)277 void calc_chi2( int vcount, int var_count, const float* vecs,
278 const float* another, Qfloat* results )
279 {
280 Mat R( 1, vcount, QFLOAT_TYPE, results );
281 double gamma = -params.gamma;
282 int j, k;
283 for( j = 0; j < vcount; j++ )
284 {
285 const float* sample = &vecs[j*var_count];
286 double chi2 = 0;
287 for(k = 0 ; k < var_count; k++ )
288 {
289 double d = sample[k]-another[k];
290 double devisor = sample[k]+another[k];
291 /// if devisor == 0, the Chi2 distance would be zero,
292 // but calculation would rise an error because of dividing by zero
293 if (devisor != 0)
294 {
295 chi2 += d*d/devisor;
296 }
297 }
298 results[j] = (Qfloat) (gamma*chi2);
299 }
300 if( vcount > 0 )
301 exp( R, R );
302 }
303
calc(int vcount,int var_count,const float * vecs,const float * another,Qfloat * results)304 void calc( int vcount, int var_count, const float* vecs,
305 const float* another, Qfloat* results ) CV_OVERRIDE
306 {
307 switch( params.kernelType )
308 {
309 case SVM::LINEAR:
310 calc_linear(vcount, var_count, vecs, another, results);
311 break;
312 case SVM::RBF:
313 calc_rbf(vcount, var_count, vecs, another, results);
314 break;
315 case SVM::POLY:
316 calc_poly(vcount, var_count, vecs, another, results);
317 break;
318 case SVM::SIGMOID:
319 calc_sigmoid(vcount, var_count, vecs, another, results);
320 break;
321 case SVM::CHI2:
322 calc_chi2(vcount, var_count, vecs, another, results);
323 break;
324 case SVM::INTER:
325 calc_intersec(vcount, var_count, vecs, another, results);
326 break;
327 default:
328 CV_Error(CV_StsBadArg, "Unknown kernel type");
329 }
330 const Qfloat max_val = (Qfloat)(FLT_MAX*1e-3);
331 for( int j = 0; j < vcount; j++ )
332 {
333 if (!(results[j] <= max_val)) // handle NaNs too
334 results[j] = max_val;
335 }
336 }
337
338 SvmParams params;
339 };
340
341
342
343 /////////////////////////////////////////////////////////////////////////
344
sortSamplesByClasses(const Mat & _samples,const Mat & _responses,vector<int> & sidx_all,vector<int> & class_ranges)345 static void sortSamplesByClasses( const Mat& _samples, const Mat& _responses,
346 vector<int>& sidx_all, vector<int>& class_ranges )
347 {
348 int i, nsamples = _samples.rows;
349 CV_Assert( _responses.isContinuous() && _responses.checkVector(1, CV_32S) == nsamples );
350
351 setRangeVector(sidx_all, nsamples);
352
353 const int* rptr = _responses.ptr<int>();
354 std::sort(sidx_all.begin(), sidx_all.end(), cmp_lt_idx<int>(rptr));
355 class_ranges.clear();
356 class_ranges.push_back(0);
357
358 for( i = 0; i < nsamples; i++ )
359 {
360 if( i == nsamples-1 || rptr[sidx_all[i]] != rptr[sidx_all[i+1]] )
361 class_ranges.push_back(i+1);
362 }
363 }
364
365 //////////////////////// SVM implementation //////////////////////////////
366
getDefaultGridPtr(int param_id)367 Ptr<ParamGrid> SVM::getDefaultGridPtr( int param_id)
368 {
369 ParamGrid grid = getDefaultGrid(param_id); // this is not a nice solution..
370 return makePtr<ParamGrid>(grid.minVal, grid.maxVal, grid.logStep);
371 }
372
getDefaultGrid(int param_id)373 ParamGrid SVM::getDefaultGrid( int param_id )
374 {
375 ParamGrid grid;
376 if( param_id == SVM::C )
377 {
378 grid.minVal = 0.1;
379 grid.maxVal = 500;
380 grid.logStep = 5; // total iterations = 5
381 }
382 else if( param_id == SVM::GAMMA )
383 {
384 grid.minVal = 1e-5;
385 grid.maxVal = 0.6;
386 grid.logStep = 15; // total iterations = 4
387 }
388 else if( param_id == SVM::P )
389 {
390 grid.minVal = 0.01;
391 grid.maxVal = 100;
392 grid.logStep = 7; // total iterations = 4
393 }
394 else if( param_id == SVM::NU )
395 {
396 grid.minVal = 0.01;
397 grid.maxVal = 0.2;
398 grid.logStep = 3; // total iterations = 3
399 }
400 else if( param_id == SVM::COEF )
401 {
402 grid.minVal = 0.1;
403 grid.maxVal = 300;
404 grid.logStep = 14; // total iterations = 3
405 }
406 else if( param_id == SVM::DEGREE )
407 {
408 grid.minVal = 0.01;
409 grid.maxVal = 4;
410 grid.logStep = 7; // total iterations = 3
411 }
412 else
413 cvError( CV_StsBadArg, "SVM::getDefaultGrid", "Invalid type of parameter "
414 "(use one of SVM::C, SVM::GAMMA et al.)", __FILE__, __LINE__ );
415 return grid;
416 }
417
418
419 class SVMImpl CV_FINAL : public SVM
420 {
421 public:
422 struct DecisionFunc
423 {
DecisionFunccv::ml::CV_FINAL::DecisionFunc424 DecisionFunc(double _rho, int _ofs) : rho(_rho), ofs(_ofs) {}
DecisionFunccv::ml::CV_FINAL::DecisionFunc425 DecisionFunc() : rho(0.), ofs(0) {}
426 double rho;
427 int ofs;
428 };
429
430 // Generalized SMO+SVMlight algorithm
431 // Solves:
432 //
433 // min [0.5(\alpha^T Q \alpha) + b^T \alpha]
434 //
435 // y^T \alpha = \delta
436 // y_i = +1 or -1
437 // 0 <= alpha_i <= Cp for y_i = 1
438 // 0 <= alpha_i <= Cn for y_i = -1
439 //
440 // Given:
441 //
442 // Q, b, y, Cp, Cn, and an initial feasible point \alpha
443 // l is the size of vectors and matrices
444 // eps is the stopping criterion
445 //
446 // solution will be put in \alpha, objective value will be put in obj
447 //
448 class Solver
449 {
450 public:
451 enum { MIN_CACHE_SIZE = (40 << 20) /* 40Mb */, MAX_CACHE_SIZE = (500 << 20) /* 500Mb */ };
452
453 typedef bool (Solver::*SelectWorkingSet)( int& i, int& j );
454 typedef Qfloat* (Solver::*GetRow)( int i, Qfloat* row, Qfloat* dst, bool existed );
455 typedef void (Solver::*CalcRho)( double& rho, double& r );
456
457 struct KernelRow
458 {
KernelRowcv::ml::CV_FINAL::Solver::KernelRow459 KernelRow() { idx = -1; prev = next = 0; }
KernelRowcv::ml::CV_FINAL::Solver::KernelRow460 KernelRow(int _idx, int _prev, int _next) : idx(_idx), prev(_prev), next(_next) {}
461 int idx;
462 int prev;
463 int next;
464 };
465
466 struct SolutionInfo
467 {
SolutionInfocv::ml::CV_FINAL::Solver::SolutionInfo468 SolutionInfo() { obj = rho = upper_bound_p = upper_bound_n = r = 0; }
469 double obj;
470 double rho;
471 double upper_bound_p;
472 double upper_bound_n;
473 double r; // for Solver_NU
474 };
475
clear()476 void clear()
477 {
478 alpha_vec = 0;
479 select_working_set_func = 0;
480 calc_rho_func = 0;
481 get_row_func = 0;
482 lru_cache.clear();
483 }
484
Solver(const Mat & _samples,const vector<schar> & _y,vector<double> & _alpha,const vector<double> & _b,double _Cp,double _Cn,const Ptr<SVM::Kernel> & _kernel,GetRow _get_row,SelectWorkingSet _select_working_set,CalcRho _calc_rho,TermCriteria _termCrit)485 Solver( const Mat& _samples, const vector<schar>& _y,
486 vector<double>& _alpha, const vector<double>& _b,
487 double _Cp, double _Cn,
488 const Ptr<SVM::Kernel>& _kernel, GetRow _get_row,
489 SelectWorkingSet _select_working_set, CalcRho _calc_rho,
490 TermCriteria _termCrit )
491 {
492 clear();
493
494 samples = _samples;
495 sample_count = samples.rows;
496 var_count = samples.cols;
497
498 y_vec = _y;
499 alpha_vec = &_alpha;
500 alpha_count = (int)alpha_vec->size();
501 b_vec = _b;
502 kernel = _kernel;
503
504 C[0] = _Cn;
505 C[1] = _Cp;
506 eps = _termCrit.epsilon;
507 max_iter = _termCrit.maxCount;
508
509 G_vec.resize(alpha_count);
510 alpha_status_vec.resize(alpha_count);
511 buf[0].resize(sample_count*2);
512 buf[1].resize(sample_count*2);
513
514 select_working_set_func = _select_working_set;
515 CV_Assert(select_working_set_func != 0);
516
517 calc_rho_func = _calc_rho;
518 CV_Assert(calc_rho_func != 0);
519
520 get_row_func = _get_row;
521 CV_Assert(get_row_func != 0);
522
523 // assume that for large training sets ~25% of Q matrix is used
524 int64 csize = (int64)sample_count*sample_count/4;
525 csize = std::max(csize, (int64)(MIN_CACHE_SIZE/sizeof(Qfloat)) );
526 csize = std::min(csize, (int64)(MAX_CACHE_SIZE/sizeof(Qfloat)) );
527 max_cache_size = (int)((csize + sample_count-1)/sample_count);
528 max_cache_size = std::min(std::max(max_cache_size, 1), sample_count);
529 cache_size = 0;
530
531 lru_cache.clear();
532 lru_cache.resize(sample_count+1, KernelRow(-1, 0, 0));
533 lru_first = lru_last = 0;
534 lru_cache_data.create(max_cache_size, sample_count, QFLOAT_TYPE);
535 }
536
get_row_base(int i,bool * _existed)537 Qfloat* get_row_base( int i, bool* _existed )
538 {
539 int i1 = i < sample_count ? i : i - sample_count;
540 KernelRow& kr = lru_cache[i1+1];
541 if( _existed )
542 *_existed = kr.idx >= 0;
543 if( kr.idx < 0 )
544 {
545 if( cache_size < max_cache_size )
546 {
547 kr.idx = cache_size;
548 cache_size++;
549 if (!lru_last)
550 lru_last = i1+1;
551 }
552 else
553 {
554 KernelRow& last = lru_cache[lru_last];
555 kr.idx = last.idx;
556 last.idx = -1;
557 lru_cache[last.prev].next = 0;
558 lru_last = last.prev;
559 last.prev = 0;
560 last.next = 0;
561 }
562 kernel->calc( sample_count, var_count, samples.ptr<float>(),
563 samples.ptr<float>(i1), lru_cache_data.ptr<Qfloat>(kr.idx) );
564 }
565 else
566 {
567 if( kr.next )
568 lru_cache[kr.next].prev = kr.prev;
569 else
570 lru_last = kr.prev;
571 if( kr.prev )
572 lru_cache[kr.prev].next = kr.next;
573 else
574 lru_first = kr.next;
575 }
576 if (lru_first)
577 lru_cache[lru_first].prev = i1+1;
578 kr.next = lru_first;
579 kr.prev = 0;
580 lru_first = i1+1;
581
582 return lru_cache_data.ptr<Qfloat>(kr.idx);
583 }
584
get_row_svc(int i,Qfloat * row,Qfloat *,bool existed)585 Qfloat* get_row_svc( int i, Qfloat* row, Qfloat*, bool existed )
586 {
587 if( !existed )
588 {
589 const schar* _y = &y_vec[0];
590 int j, len = sample_count;
591
592 if( _y[i] > 0 )
593 {
594 for( j = 0; j < len; j++ )
595 row[j] = _y[j]*row[j];
596 }
597 else
598 {
599 for( j = 0; j < len; j++ )
600 row[j] = -_y[j]*row[j];
601 }
602 }
603 return row;
604 }
605
get_row_one_class(int,Qfloat * row,Qfloat *,bool)606 Qfloat* get_row_one_class( int, Qfloat* row, Qfloat*, bool )
607 {
608 return row;
609 }
610
get_row_svr(int i,Qfloat * row,Qfloat * dst,bool)611 Qfloat* get_row_svr( int i, Qfloat* row, Qfloat* dst, bool )
612 {
613 int j, len = sample_count;
614 Qfloat* dst_pos = dst;
615 Qfloat* dst_neg = dst + len;
616 if( i >= len )
617 std::swap(dst_pos, dst_neg);
618
619 for( j = 0; j < len; j++ )
620 {
621 Qfloat t = row[j];
622 dst_pos[j] = t;
623 dst_neg[j] = -t;
624 }
625 return dst;
626 }
627
get_row(int i,float * dst)628 Qfloat* get_row( int i, float* dst )
629 {
630 bool existed = false;
631 float* row = get_row_base( i, &existed );
632 return (this->*get_row_func)( i, row, dst, existed );
633 }
634
635 #undef is_upper_bound
636 #define is_upper_bound(i) (alpha_status[i] > 0)
637
638 #undef is_lower_bound
639 #define is_lower_bound(i) (alpha_status[i] < 0)
640
641 #undef is_free
642 #define is_free(i) (alpha_status[i] == 0)
643
644 #undef get_C
645 #define get_C(i) (C[y[i]>0])
646
647 #undef update_alpha_status
648 #define update_alpha_status(i) \
649 alpha_status[i] = (schar)(alpha[i] >= get_C(i) ? 1 : alpha[i] <= 0 ? -1 : 0)
650
651 #undef reconstruct_gradient
652 #define reconstruct_gradient() /* empty for now */
653
solve_generic(SolutionInfo & si)654 bool solve_generic( SolutionInfo& si )
655 {
656 const schar* y = &y_vec[0];
657 double* alpha = &alpha_vec->at(0);
658 schar* alpha_status = &alpha_status_vec[0];
659 double* G = &G_vec[0];
660 double* b = &b_vec[0];
661
662 int iter = 0;
663 int i, j, k;
664
665 // 1. initialize gradient and alpha status
666 for( i = 0; i < alpha_count; i++ )
667 {
668 update_alpha_status(i);
669 G[i] = b[i];
670 if( fabs(G[i]) > 1e200 )
671 return false;
672 }
673
674 for( i = 0; i < alpha_count; i++ )
675 {
676 if( !is_lower_bound(i) )
677 {
678 const Qfloat *Q_i = get_row( i, &buf[0][0] );
679 double alpha_i = alpha[i];
680
681 for( j = 0; j < alpha_count; j++ )
682 G[j] += alpha_i*Q_i[j];
683 }
684 }
685
686 // 2. optimization loop
687 for(;;)
688 {
689 const Qfloat *Q_i, *Q_j;
690 double C_i, C_j;
691 double old_alpha_i, old_alpha_j, alpha_i, alpha_j;
692 double delta_alpha_i, delta_alpha_j;
693
694 #ifdef _DEBUG
695 for( i = 0; i < alpha_count; i++ )
696 {
697 if( fabs(G[i]) > 1e+300 )
698 return false;
699
700 if( fabs(alpha[i]) > 1e16 )
701 return false;
702 }
703 #endif
704
705 if( (this->*select_working_set_func)( i, j ) != 0 || iter++ >= max_iter )
706 break;
707
708 Q_i = get_row( i, &buf[0][0] );
709 Q_j = get_row( j, &buf[1][0] );
710
711 C_i = get_C(i);
712 C_j = get_C(j);
713
714 alpha_i = old_alpha_i = alpha[i];
715 alpha_j = old_alpha_j = alpha[j];
716
717 if( y[i] != y[j] )
718 {
719 double denom = Q_i[i]+Q_j[j]+2*Q_i[j];
720 double delta = (-G[i]-G[j])/MAX(fabs(denom),FLT_EPSILON);
721 double diff = alpha_i - alpha_j;
722 alpha_i += delta;
723 alpha_j += delta;
724
725 if( diff > 0 && alpha_j < 0 )
726 {
727 alpha_j = 0;
728 alpha_i = diff;
729 }
730 else if( diff <= 0 && alpha_i < 0 )
731 {
732 alpha_i = 0;
733 alpha_j = -diff;
734 }
735
736 if( diff > C_i - C_j && alpha_i > C_i )
737 {
738 alpha_i = C_i;
739 alpha_j = C_i - diff;
740 }
741 else if( diff <= C_i - C_j && alpha_j > C_j )
742 {
743 alpha_j = C_j;
744 alpha_i = C_j + diff;
745 }
746 }
747 else
748 {
749 double denom = Q_i[i]+Q_j[j]-2*Q_i[j];
750 double delta = (G[i]-G[j])/MAX(fabs(denom),FLT_EPSILON);
751 double sum = alpha_i + alpha_j;
752 alpha_i -= delta;
753 alpha_j += delta;
754
755 if( sum > C_i && alpha_i > C_i )
756 {
757 alpha_i = C_i;
758 alpha_j = sum - C_i;
759 }
760 else if( sum <= C_i && alpha_j < 0)
761 {
762 alpha_j = 0;
763 alpha_i = sum;
764 }
765
766 if( sum > C_j && alpha_j > C_j )
767 {
768 alpha_j = C_j;
769 alpha_i = sum - C_j;
770 }
771 else if( sum <= C_j && alpha_i < 0 )
772 {
773 alpha_i = 0;
774 alpha_j = sum;
775 }
776 }
777
778 // update alpha
779 alpha[i] = alpha_i;
780 alpha[j] = alpha_j;
781 update_alpha_status(i);
782 update_alpha_status(j);
783
784 // update G
785 delta_alpha_i = alpha_i - old_alpha_i;
786 delta_alpha_j = alpha_j - old_alpha_j;
787
788 for( k = 0; k < alpha_count; k++ )
789 G[k] += Q_i[k]*delta_alpha_i + Q_j[k]*delta_alpha_j;
790 }
791
792 // calculate rho
793 (this->*calc_rho_func)( si.rho, si.r );
794
795 // calculate objective value
796 for( i = 0, si.obj = 0; i < alpha_count; i++ )
797 si.obj += alpha[i] * (G[i] + b[i]);
798
799 si.obj *= 0.5;
800
801 si.upper_bound_p = C[1];
802 si.upper_bound_n = C[0];
803
804 return true;
805 }
806
807 // return 1 if already optimal, return 0 otherwise
select_working_set(int & out_i,int & out_j)808 bool select_working_set( int& out_i, int& out_j )
809 {
810 // return i,j which maximize -grad(f)^T d , under constraint
811 // if alpha_i == C, d != +1
812 // if alpha_i == 0, d != -1
813 double Gmax1 = -DBL_MAX; // max { -grad(f)_i * d | y_i*d = +1 }
814 int Gmax1_idx = -1;
815
816 double Gmax2 = -DBL_MAX; // max { -grad(f)_i * d | y_i*d = -1 }
817 int Gmax2_idx = -1;
818
819 const schar* y = &y_vec[0];
820 const schar* alpha_status = &alpha_status_vec[0];
821 const double* G = &G_vec[0];
822
823 for( int i = 0; i < alpha_count; i++ )
824 {
825 double t;
826
827 if( y[i] > 0 ) // y = +1
828 {
829 if( !is_upper_bound(i) && (t = -G[i]) > Gmax1 ) // d = +1
830 {
831 Gmax1 = t;
832 Gmax1_idx = i;
833 }
834 if( !is_lower_bound(i) && (t = G[i]) > Gmax2 ) // d = -1
835 {
836 Gmax2 = t;
837 Gmax2_idx = i;
838 }
839 }
840 else // y = -1
841 {
842 if( !is_upper_bound(i) && (t = -G[i]) > Gmax2 ) // d = +1
843 {
844 Gmax2 = t;
845 Gmax2_idx = i;
846 }
847 if( !is_lower_bound(i) && (t = G[i]) > Gmax1 ) // d = -1
848 {
849 Gmax1 = t;
850 Gmax1_idx = i;
851 }
852 }
853 }
854
855 out_i = Gmax1_idx;
856 out_j = Gmax2_idx;
857
858 return Gmax1 + Gmax2 < eps;
859 }
860
calc_rho(double & rho,double & r)861 void calc_rho( double& rho, double& r )
862 {
863 int nr_free = 0;
864 double ub = DBL_MAX, lb = -DBL_MAX, sum_free = 0;
865 const schar* y = &y_vec[0];
866 const schar* alpha_status = &alpha_status_vec[0];
867 const double* G = &G_vec[0];
868
869 for( int i = 0; i < alpha_count; i++ )
870 {
871 double yG = y[i]*G[i];
872
873 if( is_lower_bound(i) )
874 {
875 if( y[i] > 0 )
876 ub = MIN(ub,yG);
877 else
878 lb = MAX(lb,yG);
879 }
880 else if( is_upper_bound(i) )
881 {
882 if( y[i] < 0)
883 ub = MIN(ub,yG);
884 else
885 lb = MAX(lb,yG);
886 }
887 else
888 {
889 ++nr_free;
890 sum_free += yG;
891 }
892 }
893
894 rho = nr_free > 0 ? sum_free/nr_free : (ub + lb)*0.5;
895 r = 0;
896 }
897
select_working_set_nu_svm(int & out_i,int & out_j)898 bool select_working_set_nu_svm( int& out_i, int& out_j )
899 {
900 // return i,j which maximize -grad(f)^T d , under constraint
901 // if alpha_i == C, d != +1
902 // if alpha_i == 0, d != -1
903 double Gmax1 = -DBL_MAX; // max { -grad(f)_i * d | y_i = +1, d = +1 }
904 int Gmax1_idx = -1;
905
906 double Gmax2 = -DBL_MAX; // max { -grad(f)_i * d | y_i = +1, d = -1 }
907 int Gmax2_idx = -1;
908
909 double Gmax3 = -DBL_MAX; // max { -grad(f)_i * d | y_i = -1, d = +1 }
910 int Gmax3_idx = -1;
911
912 double Gmax4 = -DBL_MAX; // max { -grad(f)_i * d | y_i = -1, d = -1 }
913 int Gmax4_idx = -1;
914
915 const schar* y = &y_vec[0];
916 const schar* alpha_status = &alpha_status_vec[0];
917 const double* G = &G_vec[0];
918
919 for( int i = 0; i < alpha_count; i++ )
920 {
921 double t;
922
923 if( y[i] > 0 ) // y == +1
924 {
925 if( !is_upper_bound(i) && (t = -G[i]) > Gmax1 ) // d = +1
926 {
927 Gmax1 = t;
928 Gmax1_idx = i;
929 }
930 if( !is_lower_bound(i) && (t = G[i]) > Gmax2 ) // d = -1
931 {
932 Gmax2 = t;
933 Gmax2_idx = i;
934 }
935 }
936 else // y == -1
937 {
938 if( !is_upper_bound(i) && (t = -G[i]) > Gmax3 ) // d = +1
939 {
940 Gmax3 = t;
941 Gmax3_idx = i;
942 }
943 if( !is_lower_bound(i) && (t = G[i]) > Gmax4 ) // d = -1
944 {
945 Gmax4 = t;
946 Gmax4_idx = i;
947 }
948 }
949 }
950
951 if( MAX(Gmax1 + Gmax2, Gmax3 + Gmax4) < eps )
952 return 1;
953
954 if( Gmax1 + Gmax2 > Gmax3 + Gmax4 )
955 {
956 out_i = Gmax1_idx;
957 out_j = Gmax2_idx;
958 }
959 else
960 {
961 out_i = Gmax3_idx;
962 out_j = Gmax4_idx;
963 }
964 return 0;
965 }
966
calc_rho_nu_svm(double & rho,double & r)967 void calc_rho_nu_svm( double& rho, double& r )
968 {
969 int nr_free1 = 0, nr_free2 = 0;
970 double ub1 = DBL_MAX, ub2 = DBL_MAX;
971 double lb1 = -DBL_MAX, lb2 = -DBL_MAX;
972 double sum_free1 = 0, sum_free2 = 0;
973
974 const schar* y = &y_vec[0];
975 const schar* alpha_status = &alpha_status_vec[0];
976 const double* G = &G_vec[0];
977
978 for( int i = 0; i < alpha_count; i++ )
979 {
980 double G_i = G[i];
981 if( y[i] > 0 )
982 {
983 if( is_lower_bound(i) )
984 ub1 = MIN( ub1, G_i );
985 else if( is_upper_bound(i) )
986 lb1 = MAX( lb1, G_i );
987 else
988 {
989 ++nr_free1;
990 sum_free1 += G_i;
991 }
992 }
993 else
994 {
995 if( is_lower_bound(i) )
996 ub2 = MIN( ub2, G_i );
997 else if( is_upper_bound(i) )
998 lb2 = MAX( lb2, G_i );
999 else
1000 {
1001 ++nr_free2;
1002 sum_free2 += G_i;
1003 }
1004 }
1005 }
1006
1007 double r1 = nr_free1 > 0 ? sum_free1/nr_free1 : (ub1 + lb1)*0.5;
1008 double r2 = nr_free2 > 0 ? sum_free2/nr_free2 : (ub2 + lb2)*0.5;
1009
1010 rho = (r1 - r2)*0.5;
1011 r = (r1 + r2)*0.5;
1012 }
1013
1014 /*
1015 ///////////////////////// construct and solve various formulations ///////////////////////
1016 */
solve_c_svc(const Mat & _samples,const vector<schar> & _y,double _Cp,double _Cn,const Ptr<SVM::Kernel> & _kernel,vector<double> & _alpha,SolutionInfo & _si,TermCriteria termCrit)1017 static bool solve_c_svc( const Mat& _samples, const vector<schar>& _y,
1018 double _Cp, double _Cn, const Ptr<SVM::Kernel>& _kernel,
1019 vector<double>& _alpha, SolutionInfo& _si, TermCriteria termCrit )
1020 {
1021 int sample_count = _samples.rows;
1022
1023 _alpha.assign(sample_count, 0.);
1024 vector<double> _b(sample_count, -1.);
1025
1026 Solver solver( _samples, _y, _alpha, _b, _Cp, _Cn, _kernel,
1027 &Solver::get_row_svc,
1028 &Solver::select_working_set,
1029 &Solver::calc_rho,
1030 termCrit );
1031
1032 if( !solver.solve_generic( _si ))
1033 return false;
1034
1035 for( int i = 0; i < sample_count; i++ )
1036 _alpha[i] *= _y[i];
1037
1038 return true;
1039 }
1040
1041
solve_nu_svc(const Mat & _samples,const vector<schar> & _y,double nu,const Ptr<SVM::Kernel> & _kernel,vector<double> & _alpha,SolutionInfo & _si,TermCriteria termCrit)1042 static bool solve_nu_svc( const Mat& _samples, const vector<schar>& _y,
1043 double nu, const Ptr<SVM::Kernel>& _kernel,
1044 vector<double>& _alpha, SolutionInfo& _si,
1045 TermCriteria termCrit )
1046 {
1047 int sample_count = _samples.rows;
1048
1049 _alpha.resize(sample_count);
1050 vector<double> _b(sample_count, 0.);
1051
1052 double sum_pos = nu * sample_count * 0.5;
1053 double sum_neg = nu * sample_count * 0.5;
1054
1055 for( int i = 0; i < sample_count; i++ )
1056 {
1057 double a;
1058 if( _y[i] > 0 )
1059 {
1060 a = std::min(1.0, sum_pos);
1061 sum_pos -= a;
1062 }
1063 else
1064 {
1065 a = std::min(1.0, sum_neg);
1066 sum_neg -= a;
1067 }
1068 _alpha[i] = a;
1069 }
1070
1071 Solver solver( _samples, _y, _alpha, _b, 1., 1., _kernel,
1072 &Solver::get_row_svc,
1073 &Solver::select_working_set_nu_svm,
1074 &Solver::calc_rho_nu_svm,
1075 termCrit );
1076
1077 if( !solver.solve_generic( _si ))
1078 return false;
1079
1080 double inv_r = 1./_si.r;
1081
1082 for( int i = 0; i < sample_count; i++ )
1083 _alpha[i] *= _y[i]*inv_r;
1084
1085 _si.rho *= inv_r;
1086 _si.obj *= (inv_r*inv_r);
1087 _si.upper_bound_p = inv_r;
1088 _si.upper_bound_n = inv_r;
1089
1090 return true;
1091 }
1092
solve_one_class(const Mat & _samples,double nu,const Ptr<SVM::Kernel> & _kernel,vector<double> & _alpha,SolutionInfo & _si,TermCriteria termCrit)1093 static bool solve_one_class( const Mat& _samples, double nu,
1094 const Ptr<SVM::Kernel>& _kernel,
1095 vector<double>& _alpha, SolutionInfo& _si,
1096 TermCriteria termCrit )
1097 {
1098 int sample_count = _samples.rows;
1099 vector<schar> _y(sample_count, 1);
1100 vector<double> _b(sample_count, 0.);
1101
1102 int i, n = cvRound( nu*sample_count );
1103
1104 _alpha.resize(sample_count);
1105 for( i = 0; i < sample_count; i++ )
1106 _alpha[i] = i < n ? 1 : 0;
1107
1108 if( n < sample_count )
1109 _alpha[n] = nu * sample_count - n;
1110 else
1111 _alpha[n-1] = nu * sample_count - (n-1);
1112
1113 Solver solver( _samples, _y, _alpha, _b, 1., 1., _kernel,
1114 &Solver::get_row_one_class,
1115 &Solver::select_working_set,
1116 &Solver::calc_rho,
1117 termCrit );
1118
1119 return solver.solve_generic(_si);
1120 }
1121
solve_eps_svr(const Mat & _samples,const vector<float> & _yf,double p,double C,const Ptr<SVM::Kernel> & _kernel,vector<double> & _alpha,SolutionInfo & _si,TermCriteria termCrit)1122 static bool solve_eps_svr( const Mat& _samples, const vector<float>& _yf,
1123 double p, double C, const Ptr<SVM::Kernel>& _kernel,
1124 vector<double>& _alpha, SolutionInfo& _si,
1125 TermCriteria termCrit )
1126 {
1127 int sample_count = _samples.rows;
1128 int alpha_count = sample_count*2;
1129
1130 CV_Assert( (int)_yf.size() == sample_count );
1131
1132 _alpha.assign(alpha_count, 0.);
1133 vector<schar> _y(alpha_count);
1134 vector<double> _b(alpha_count);
1135
1136 for( int i = 0; i < sample_count; i++ )
1137 {
1138 _b[i] = p - _yf[i];
1139 _y[i] = 1;
1140
1141 _b[i+sample_count] = p + _yf[i];
1142 _y[i+sample_count] = -1;
1143 }
1144
1145 Solver solver( _samples, _y, _alpha, _b, C, C, _kernel,
1146 &Solver::get_row_svr,
1147 &Solver::select_working_set,
1148 &Solver::calc_rho,
1149 termCrit );
1150
1151 if( !solver.solve_generic( _si ))
1152 return false;
1153
1154 for( int i = 0; i < sample_count; i++ )
1155 _alpha[i] -= _alpha[i+sample_count];
1156
1157 return true;
1158 }
1159
1160
solve_nu_svr(const Mat & _samples,const vector<float> & _yf,double nu,double C,const Ptr<SVM::Kernel> & _kernel,vector<double> & _alpha,SolutionInfo & _si,TermCriteria termCrit)1161 static bool solve_nu_svr( const Mat& _samples, const vector<float>& _yf,
1162 double nu, double C, const Ptr<SVM::Kernel>& _kernel,
1163 vector<double>& _alpha, SolutionInfo& _si,
1164 TermCriteria termCrit )
1165 {
1166 int sample_count = _samples.rows;
1167 int alpha_count = sample_count*2;
1168 double sum = C * nu * sample_count * 0.5;
1169
1170 CV_Assert( (int)_yf.size() == sample_count );
1171
1172 _alpha.resize(alpha_count);
1173 vector<schar> _y(alpha_count);
1174 vector<double> _b(alpha_count);
1175
1176 for( int i = 0; i < sample_count; i++ )
1177 {
1178 _alpha[i] = _alpha[i + sample_count] = std::min(sum, C);
1179 sum -= _alpha[i];
1180
1181 _b[i] = -_yf[i];
1182 _y[i] = 1;
1183
1184 _b[i + sample_count] = _yf[i];
1185 _y[i + sample_count] = -1;
1186 }
1187
1188 Solver solver( _samples, _y, _alpha, _b, 1., 1., _kernel,
1189 &Solver::get_row_svr,
1190 &Solver::select_working_set_nu_svm,
1191 &Solver::calc_rho_nu_svm,
1192 termCrit );
1193
1194 if( !solver.solve_generic( _si ))
1195 return false;
1196
1197 for( int i = 0; i < sample_count; i++ )
1198 _alpha[i] -= _alpha[i+sample_count];
1199
1200 return true;
1201 }
1202
1203 int sample_count;
1204 int var_count;
1205 int cache_size;
1206 int max_cache_size;
1207 Mat samples;
1208 SvmParams params;
1209 vector<KernelRow> lru_cache;
1210 int lru_first;
1211 int lru_last;
1212 Mat lru_cache_data;
1213
1214 int alpha_count;
1215
1216 vector<double> G_vec;
1217 vector<double>* alpha_vec;
1218 vector<schar> y_vec;
1219 // -1 - lower bound, 0 - free, 1 - upper bound
1220 vector<schar> alpha_status_vec;
1221 vector<double> b_vec;
1222
1223 vector<Qfloat> buf[2];
1224 double eps;
1225 int max_iter;
1226 double C[2]; // C[0] == Cn, C[1] == Cp
1227 Ptr<SVM::Kernel> kernel;
1228
1229 SelectWorkingSet select_working_set_func;
1230 CalcRho calc_rho_func;
1231 GetRow get_row_func;
1232 };
1233
1234 //////////////////////////////////////////////////////////////////////////////////////////
SVMImpl()1235 SVMImpl()
1236 {
1237 clear();
1238 checkParams();
1239 }
1240
~SVMImpl()1241 ~SVMImpl()
1242 {
1243 clear();
1244 }
1245
clear()1246 void clear() CV_OVERRIDE
1247 {
1248 decision_func.clear();
1249 df_alpha.clear();
1250 df_index.clear();
1251 sv.release();
1252 uncompressed_sv.release();
1253 }
1254
getUncompressedSupportVectors() const1255 Mat getUncompressedSupportVectors() const CV_OVERRIDE
1256 {
1257 return uncompressed_sv;
1258 }
1259
getSupportVectors() const1260 Mat getSupportVectors() const CV_OVERRIDE
1261 {
1262 return sv;
1263 }
1264
getType() const1265 inline int getType() const CV_OVERRIDE { return params.svmType; }
setType(int val)1266 inline void setType(int val) CV_OVERRIDE { params.svmType = val; }
getGamma() const1267 inline double getGamma() const CV_OVERRIDE { return params.gamma; }
setGamma(double val)1268 inline void setGamma(double val) CV_OVERRIDE { params.gamma = val; }
getCoef0() const1269 inline double getCoef0() const CV_OVERRIDE { return params.coef0; }
setCoef0(double val)1270 inline void setCoef0(double val) CV_OVERRIDE { params.coef0 = val; }
getDegree() const1271 inline double getDegree() const CV_OVERRIDE { return params.degree; }
setDegree(double val)1272 inline void setDegree(double val) CV_OVERRIDE { params.degree = val; }
getC() const1273 inline double getC() const CV_OVERRIDE { return params.C; }
setC(double val)1274 inline void setC(double val) CV_OVERRIDE { params.C = val; }
getNu() const1275 inline double getNu() const CV_OVERRIDE { return params.nu; }
setNu(double val)1276 inline void setNu(double val) CV_OVERRIDE { params.nu = val; }
getP() const1277 inline double getP() const CV_OVERRIDE { return params.p; }
setP(double val)1278 inline void setP(double val) CV_OVERRIDE { params.p = val; }
getClassWeights() const1279 inline cv::Mat getClassWeights() const CV_OVERRIDE { return params.classWeights; }
setClassWeights(const cv::Mat & val)1280 inline void setClassWeights(const cv::Mat& val) CV_OVERRIDE { params.classWeights = val; }
getTermCriteria() const1281 inline cv::TermCriteria getTermCriteria() const CV_OVERRIDE { return params.termCrit; }
setTermCriteria(const cv::TermCriteria & val)1282 inline void setTermCriteria(const cv::TermCriteria& val) CV_OVERRIDE { params.termCrit = val; }
1283
getKernelType() const1284 int getKernelType() const CV_OVERRIDE { return params.kernelType; }
setKernel(int kernelType)1285 void setKernel(int kernelType) CV_OVERRIDE
1286 {
1287 params.kernelType = kernelType;
1288 if (kernelType != CUSTOM)
1289 kernel = makePtr<SVMKernelImpl>(params);
1290 }
1291
setCustomKernel(const Ptr<Kernel> & _kernel)1292 void setCustomKernel(const Ptr<Kernel> &_kernel) CV_OVERRIDE
1293 {
1294 params.kernelType = CUSTOM;
1295 kernel = _kernel;
1296 }
1297
checkParams()1298 void checkParams()
1299 {
1300 int kernelType = params.kernelType;
1301 if (kernelType != CUSTOM)
1302 {
1303 if( kernelType != LINEAR && kernelType != POLY &&
1304 kernelType != SIGMOID && kernelType != RBF &&
1305 kernelType != INTER && kernelType != CHI2)
1306 CV_Error( CV_StsBadArg, "Unknown/unsupported kernel type" );
1307
1308 if( kernelType == LINEAR )
1309 params.gamma = 1;
1310 else if( params.gamma <= 0 )
1311 CV_Error( CV_StsOutOfRange, "gamma parameter of the kernel must be positive" );
1312
1313 if( kernelType != SIGMOID && kernelType != POLY )
1314 params.coef0 = 0;
1315
1316 if( kernelType != POLY )
1317 params.degree = 0;
1318 else if( params.degree <= 0 )
1319 CV_Error( CV_StsOutOfRange, "The kernel parameter <degree> must be positive" );
1320
1321 kernel = makePtr<SVMKernelImpl>(params);
1322 }
1323 else
1324 {
1325 if (!kernel)
1326 CV_Error( CV_StsBadArg, "Custom kernel is not set" );
1327 }
1328
1329 int svmType = params.svmType;
1330
1331 if( svmType != C_SVC && svmType != NU_SVC &&
1332 svmType != ONE_CLASS && svmType != EPS_SVR &&
1333 svmType != NU_SVR )
1334 CV_Error( CV_StsBadArg, "Unknown/unsupported SVM type" );
1335
1336 if( svmType == ONE_CLASS || svmType == NU_SVC )
1337 params.C = 0;
1338 else if( params.C <= 0 )
1339 CV_Error( CV_StsOutOfRange, "The parameter C must be positive" );
1340
1341 if( svmType == C_SVC || svmType == EPS_SVR )
1342 params.nu = 0;
1343 else if( params.nu <= 0 || params.nu >= 1 )
1344 CV_Error( CV_StsOutOfRange, "The parameter nu must be between 0 and 1" );
1345
1346 if( svmType != EPS_SVR )
1347 params.p = 0;
1348 else if( params.p <= 0 )
1349 CV_Error( CV_StsOutOfRange, "The parameter p must be positive" );
1350
1351 if( svmType != C_SVC )
1352 params.classWeights.release();
1353
1354 if( !(params.termCrit.type & TermCriteria::EPS) )
1355 params.termCrit.epsilon = DBL_EPSILON;
1356 params.termCrit.epsilon = std::max(params.termCrit.epsilon, DBL_EPSILON);
1357 if( !(params.termCrit.type & TermCriteria::COUNT) )
1358 params.termCrit.maxCount = INT_MAX;
1359 params.termCrit.maxCount = std::max(params.termCrit.maxCount, 1);
1360 }
1361
setParams(const SvmParams & _params)1362 void setParams( const SvmParams& _params)
1363 {
1364 params = _params;
1365 checkParams();
1366 }
1367
getSVCount(int i) const1368 int getSVCount(int i) const
1369 {
1370 return (i < (int)(decision_func.size()-1) ? decision_func[i+1].ofs :
1371 (int)df_index.size()) - decision_func[i].ofs;
1372 }
1373
do_train(const Mat & _samples,const Mat & _responses)1374 bool do_train( const Mat& _samples, const Mat& _responses )
1375 {
1376 int svmType = params.svmType;
1377 int i, j, k, sample_count = _samples.rows;
1378 vector<double> _alpha;
1379 Solver::SolutionInfo sinfo;
1380
1381 CV_Assert( _samples.type() == CV_32F );
1382 var_count = _samples.cols;
1383
1384 if( svmType == ONE_CLASS || svmType == EPS_SVR || svmType == NU_SVR )
1385 {
1386 int sv_count = 0;
1387 decision_func.clear();
1388
1389 vector<float> _yf;
1390 if( !_responses.empty() )
1391 _responses.convertTo(_yf, CV_32F);
1392
1393 bool ok =
1394 svmType == ONE_CLASS ? Solver::solve_one_class( _samples, params.nu, kernel, _alpha, sinfo, params.termCrit ) :
1395 svmType == EPS_SVR ? Solver::solve_eps_svr( _samples, _yf, params.p, params.C, kernel, _alpha, sinfo, params.termCrit ) :
1396 svmType == NU_SVR ? Solver::solve_nu_svr( _samples, _yf, params.nu, params.C, kernel, _alpha, sinfo, params.termCrit ) : false;
1397
1398 if( !ok )
1399 return false;
1400
1401 for( i = 0; i < sample_count; i++ )
1402 sv_count += fabs(_alpha[i]) > 0;
1403
1404 CV_Assert(sv_count != 0);
1405
1406 sv.create(sv_count, _samples.cols, CV_32F);
1407 df_alpha.resize(sv_count);
1408 df_index.resize(sv_count);
1409
1410 for( i = k = 0; i < sample_count; i++ )
1411 {
1412 if( std::abs(_alpha[i]) > 0 )
1413 {
1414 _samples.row(i).copyTo(sv.row(k));
1415 df_alpha[k] = _alpha[i];
1416 df_index[k] = k;
1417 k++;
1418 }
1419 }
1420
1421 decision_func.push_back(DecisionFunc(sinfo.rho, 0));
1422 }
1423 else
1424 {
1425 int class_count = (int)class_labels.total();
1426 vector<int> svidx, sidx, sidx_all, sv_tab(sample_count, 0);
1427 Mat temp_samples, class_weights;
1428 vector<int> class_ranges;
1429 vector<schar> temp_y;
1430 double nu = params.nu;
1431 CV_Assert( svmType == C_SVC || svmType == NU_SVC );
1432
1433 if( svmType == C_SVC && !params.classWeights.empty() )
1434 {
1435 const Mat cw = params.classWeights;
1436
1437 if( (cw.cols != 1 && cw.rows != 1) ||
1438 (int)cw.total() != class_count ||
1439 (cw.type() != CV_32F && cw.type() != CV_64F) )
1440 CV_Error( CV_StsBadArg, "params.class_weights must be 1d floating-point vector "
1441 "containing as many elements as the number of classes" );
1442
1443 cw.convertTo(class_weights, CV_64F, params.C);
1444 //normalize(cw, class_weights, params.C, 0, NORM_L1, CV_64F);
1445 }
1446
1447 decision_func.clear();
1448 df_alpha.clear();
1449 df_index.clear();
1450
1451 sortSamplesByClasses( _samples, _responses, sidx_all, class_ranges );
1452
1453 //check that while cross-validation there were the samples from all the classes
1454 if ((int)class_ranges.size() < class_count + 1)
1455 CV_Error( CV_StsBadArg, "While cross-validation one or more of the classes have "
1456 "been fell out of the sample. Try to reduce <Params::k_fold>" );
1457
1458 if( svmType == NU_SVC )
1459 {
1460 // check if nu is feasible
1461 for( i = 0; i < class_count; i++ )
1462 {
1463 int ci = class_ranges[i+1] - class_ranges[i];
1464 for( j = i+1; j< class_count; j++ )
1465 {
1466 int cj = class_ranges[j+1] - class_ranges[j];
1467 if( nu*(ci + cj)*0.5 > std::min( ci, cj ) )
1468 // TODO: add some diagnostic
1469 return false;
1470 }
1471 }
1472 }
1473
1474 size_t samplesize = _samples.cols*_samples.elemSize();
1475
1476 // train n*(n-1)/2 classifiers
1477 for( i = 0; i < class_count; i++ )
1478 {
1479 for( j = i+1; j < class_count; j++ )
1480 {
1481 int si = class_ranges[i], ci = class_ranges[i+1] - si;
1482 int sj = class_ranges[j], cj = class_ranges[j+1] - sj;
1483 double Cp = params.C, Cn = Cp;
1484
1485 temp_samples.create(ci + cj, _samples.cols, _samples.type());
1486 sidx.resize(ci + cj);
1487 temp_y.resize(ci + cj);
1488
1489 // form input for the binary classification problem
1490 for( k = 0; k < ci+cj; k++ )
1491 {
1492 int idx = k < ci ? si+k : sj+k-ci;
1493 memcpy(temp_samples.ptr(k), _samples.ptr(sidx_all[idx]), samplesize);
1494 sidx[k] = sidx_all[idx];
1495 temp_y[k] = k < ci ? 1 : -1;
1496 }
1497
1498 if( !class_weights.empty() )
1499 {
1500 Cp = class_weights.at<double>(i);
1501 Cn = class_weights.at<double>(j);
1502 }
1503
1504 DecisionFunc df;
1505 bool ok = params.svmType == C_SVC ?
1506 Solver::solve_c_svc( temp_samples, temp_y, Cp, Cn,
1507 kernel, _alpha, sinfo, params.termCrit ) :
1508 params.svmType == NU_SVC ?
1509 Solver::solve_nu_svc( temp_samples, temp_y, params.nu,
1510 kernel, _alpha, sinfo, params.termCrit ) :
1511 false;
1512 if( !ok )
1513 return false;
1514 df.rho = sinfo.rho;
1515 df.ofs = (int)df_index.size();
1516 decision_func.push_back(df);
1517
1518 for( k = 0; k < ci + cj; k++ )
1519 {
1520 if( std::abs(_alpha[k]) > 0 )
1521 {
1522 int idx = k < ci ? si+k : sj+k-ci;
1523 sv_tab[sidx_all[idx]] = 1;
1524 df_index.push_back(sidx_all[idx]);
1525 df_alpha.push_back(_alpha[k]);
1526 }
1527 }
1528 }
1529 }
1530
1531 // allocate support vectors and initialize sv_tab
1532 for( i = 0, k = 0; i < sample_count; i++ )
1533 {
1534 if( sv_tab[i] )
1535 sv_tab[i] = ++k;
1536 }
1537
1538 int sv_total = k;
1539 sv.create(sv_total, _samples.cols, _samples.type());
1540
1541 for( i = 0; i < sample_count; i++ )
1542 {
1543 if( !sv_tab[i] )
1544 continue;
1545 memcpy(sv.ptr(sv_tab[i]-1), _samples.ptr(i), samplesize);
1546 }
1547
1548 // set sv pointers
1549 int n = (int)df_index.size();
1550 for( i = 0; i < n; i++ )
1551 {
1552 CV_Assert( sv_tab[df_index[i]] > 0 );
1553 df_index[i] = sv_tab[df_index[i]] - 1;
1554 }
1555 }
1556
1557 optimize_linear_svm();
1558
1559 return true;
1560 }
1561
optimize_linear_svm()1562 void optimize_linear_svm()
1563 {
1564 // we optimize only linear SVM: compress all the support vectors into one.
1565 if( params.kernelType != LINEAR )
1566 return;
1567
1568 int i, df_count = (int)decision_func.size();
1569
1570 for( i = 0; i < df_count; i++ )
1571 {
1572 if( getSVCount(i) != 1 )
1573 break;
1574 }
1575
1576 // if every decision functions uses a single support vector;
1577 // it's already compressed. skip it then.
1578 if( i == df_count )
1579 return;
1580
1581 AutoBuffer<double> vbuf(var_count);
1582 double* v = vbuf.data();
1583 Mat new_sv(df_count, var_count, CV_32F);
1584
1585 vector<DecisionFunc> new_df;
1586
1587 for( i = 0; i < df_count; i++ )
1588 {
1589 float* dst = new_sv.ptr<float>(i);
1590 memset(v, 0, var_count*sizeof(v[0]));
1591 int j, k, sv_count = getSVCount(i);
1592 const DecisionFunc& df = decision_func[i];
1593 const int* sv_index = &df_index[df.ofs];
1594 const double* sv_alpha = &df_alpha[df.ofs];
1595 for( j = 0; j < sv_count; j++ )
1596 {
1597 const float* src = sv.ptr<float>(sv_index[j]);
1598 double a = sv_alpha[j];
1599 for( k = 0; k < var_count; k++ )
1600 v[k] += src[k]*a;
1601 }
1602 for( k = 0; k < var_count; k++ )
1603 dst[k] = (float)v[k];
1604 new_df.push_back(DecisionFunc(df.rho, i));
1605 }
1606
1607 setRangeVector(df_index, df_count);
1608 df_alpha.assign(df_count, 1.);
1609 sv.copyTo(uncompressed_sv);
1610 std::swap(sv, new_sv);
1611 std::swap(decision_func, new_df);
1612 }
1613
train(const Ptr<TrainData> & data,int)1614 bool train( const Ptr<TrainData>& data, int ) CV_OVERRIDE
1615 {
1616 CV_Assert(!data.empty());
1617 clear();
1618
1619 checkParams();
1620
1621 int svmType = params.svmType;
1622 Mat samples = data->getTrainSamples();
1623 Mat responses;
1624
1625 if( svmType == C_SVC || svmType == NU_SVC )
1626 {
1627 responses = data->getTrainNormCatResponses();
1628 if( responses.empty() )
1629 CV_Error(CV_StsBadArg, "in the case of classification problem the responses must be categorical; "
1630 "either specify varType when creating TrainData, or pass integer responses");
1631 class_labels = data->getClassLabels();
1632 }
1633 else
1634 responses = data->getTrainResponses();
1635
1636 if( !do_train( samples, responses ))
1637 {
1638 clear();
1639 return false;
1640 }
1641
1642 return true;
1643 }
1644
1645 class TrainAutoBody : public ParallelLoopBody
1646 {
1647 public:
TrainAutoBody(const vector<SvmParams> & _parameters,const cv::Mat & _samples,const cv::Mat & _responses,const cv::Mat & _labels,const vector<int> & _sidx,bool _is_classification,int _k_fold,std::vector<double> & _result)1648 TrainAutoBody(const vector<SvmParams>& _parameters,
1649 const cv::Mat& _samples,
1650 const cv::Mat& _responses,
1651 const cv::Mat& _labels,
1652 const vector<int>& _sidx,
1653 bool _is_classification,
1654 int _k_fold,
1655 std::vector<double>& _result) :
1656 parameters(_parameters), samples(_samples), responses(_responses), labels(_labels),
1657 sidx(_sidx), is_classification(_is_classification), k_fold(_k_fold), result(_result)
1658 {}
1659
operator ()(const cv::Range & range) const1660 void operator()( const cv::Range& range ) const CV_OVERRIDE
1661 {
1662 int sample_count = samples.rows;
1663 int var_count_ = samples.cols;
1664 size_t sample_size = var_count_*samples.elemSize();
1665
1666 int test_sample_count = (sample_count + k_fold/2)/k_fold;
1667 int train_sample_count = sample_count - test_sample_count;
1668
1669 // Use a local instance
1670 cv::Ptr<SVMImpl> svm = makePtr<SVMImpl>();
1671 svm->class_labels = labels;
1672
1673 int rtype = responses.type();
1674
1675 Mat temp_train_samples(train_sample_count, var_count_, CV_32F);
1676 Mat temp_test_samples(test_sample_count, var_count_, CV_32F);
1677 Mat temp_train_responses(train_sample_count, 1, rtype);
1678 Mat temp_test_responses;
1679
1680 for( int p = range.start; p < range.end; p++ )
1681 {
1682 svm->setParams(parameters[p]);
1683
1684 double error = 0;
1685 for( int k = 0; k < k_fold; k++ )
1686 {
1687 int start = (k*sample_count + k_fold/2)/k_fold;
1688 for( int i = 0; i < train_sample_count; i++ )
1689 {
1690 int j = sidx[(i+start)%sample_count];
1691 memcpy(temp_train_samples.ptr(i), samples.ptr(j), sample_size);
1692 if( is_classification )
1693 temp_train_responses.at<int>(i) = responses.at<int>(j);
1694 else if( !responses.empty() )
1695 temp_train_responses.at<float>(i) = responses.at<float>(j);
1696 }
1697
1698 // Train SVM on <train_size> samples
1699 if( !svm->do_train( temp_train_samples, temp_train_responses ))
1700 continue;
1701
1702 for( int i = 0; i < test_sample_count; i++ )
1703 {
1704 int j = sidx[(i+start+train_sample_count) % sample_count];
1705 memcpy(temp_test_samples.ptr(i), samples.ptr(j), sample_size);
1706 }
1707
1708 svm->predict(temp_test_samples, temp_test_responses, 0);
1709 for( int i = 0; i < test_sample_count; i++ )
1710 {
1711 float val = temp_test_responses.at<float>(i);
1712 int j = sidx[(i+start+train_sample_count) % sample_count];
1713 if( is_classification )
1714 error += (float)(val != responses.at<int>(j));
1715 else
1716 {
1717 val -= responses.at<float>(j);
1718 error += val*val;
1719 }
1720 }
1721 }
1722
1723 result[p] = error;
1724 }
1725 }
1726
1727 private:
1728 const vector<SvmParams>& parameters;
1729 const cv::Mat& samples;
1730 const cv::Mat& responses;
1731 const cv::Mat& labels;
1732 const vector<int>& sidx;
1733 bool is_classification;
1734 int k_fold;
1735 std::vector<double>& result;
1736 };
1737
trainAuto(const Ptr<TrainData> & data,int k_fold,ParamGrid C_grid,ParamGrid gamma_grid,ParamGrid p_grid,ParamGrid nu_grid,ParamGrid coef_grid,ParamGrid degree_grid,bool balanced)1738 bool trainAuto( const Ptr<TrainData>& data, int k_fold,
1739 ParamGrid C_grid, ParamGrid gamma_grid, ParamGrid p_grid,
1740 ParamGrid nu_grid, ParamGrid coef_grid, ParamGrid degree_grid,
1741 bool balanced ) CV_OVERRIDE
1742 {
1743 CV_Assert(!data.empty());
1744 checkParams();
1745
1746 int svmType = params.svmType;
1747 RNG rng((uint64)-1);
1748
1749 if( svmType == ONE_CLASS )
1750 // current implementation of "auto" svm does not support the 1-class case.
1751 return train( data, 0 );
1752
1753 clear();
1754
1755 CV_Assert( k_fold >= 2 );
1756
1757 // All the parameters except, possibly, <coef0> are positive.
1758 // <coef0> is nonnegative
1759 #define CHECK_GRID(grid, param) \
1760 if( grid.logStep <= 1 ) \
1761 { \
1762 grid.minVal = grid.maxVal = params.param; \
1763 grid.logStep = 10; \
1764 } \
1765 else \
1766 checkParamGrid(grid)
1767
1768 CHECK_GRID(C_grid, C);
1769 CHECK_GRID(gamma_grid, gamma);
1770 CHECK_GRID(p_grid, p);
1771 CHECK_GRID(nu_grid, nu);
1772 CHECK_GRID(coef_grid, coef0);
1773 CHECK_GRID(degree_grid, degree);
1774
1775 // these parameters are not used:
1776 if( params.kernelType != POLY )
1777 degree_grid.minVal = degree_grid.maxVal = params.degree;
1778 if( params.kernelType == LINEAR )
1779 gamma_grid.minVal = gamma_grid.maxVal = params.gamma;
1780 if( params.kernelType != POLY && params.kernelType != SIGMOID )
1781 coef_grid.minVal = coef_grid.maxVal = params.coef0;
1782 if( svmType == NU_SVC || svmType == ONE_CLASS )
1783 C_grid.minVal = C_grid.maxVal = params.C;
1784 if( svmType == C_SVC || svmType == EPS_SVR )
1785 nu_grid.minVal = nu_grid.maxVal = params.nu;
1786 if( svmType != EPS_SVR )
1787 p_grid.minVal = p_grid.maxVal = params.p;
1788
1789 Mat samples = data->getTrainSamples();
1790 Mat responses;
1791 bool is_classification = false;
1792 Mat class_labels0;
1793 int class_count = (int)class_labels.total();
1794
1795 if( svmType == C_SVC || svmType == NU_SVC )
1796 {
1797 responses = data->getTrainNormCatResponses();
1798 class_labels = data->getClassLabels();
1799 class_count = (int)class_labels.total();
1800 is_classification = true;
1801
1802 vector<int> temp_class_labels;
1803 setRangeVector(temp_class_labels, class_count);
1804
1805 // temporarily replace class labels with 0, 1, ..., NCLASSES-1
1806 class_labels0 = class_labels;
1807 class_labels = Mat(temp_class_labels).clone();
1808 }
1809 else
1810 responses = data->getTrainResponses();
1811
1812 CV_Assert(samples.type() == CV_32F);
1813
1814 int sample_count = samples.rows;
1815 var_count = samples.cols;
1816
1817 vector<int> sidx;
1818 setRangeVector(sidx, sample_count);
1819
1820 // randomly permute training samples
1821 for( int i = 0; i < sample_count; i++ )
1822 {
1823 int i1 = rng.uniform(0, sample_count);
1824 int i2 = rng.uniform(0, sample_count);
1825 std::swap(sidx[i1], sidx[i2]);
1826 }
1827
1828 if( is_classification && class_count == 2 && balanced )
1829 {
1830 // reshuffle the training set in such a way that
1831 // instances of each class are divided more or less evenly
1832 // between the k_fold parts.
1833 vector<int> sidx0, sidx1;
1834
1835 for( int i = 0; i < sample_count; i++ )
1836 {
1837 if( responses.at<int>(sidx[i]) == 0 )
1838 sidx0.push_back(sidx[i]);
1839 else
1840 sidx1.push_back(sidx[i]);
1841 }
1842
1843 int n0 = (int)sidx0.size(), n1 = (int)sidx1.size();
1844 int a0 = 0, a1 = 0;
1845 sidx.clear();
1846 for( int k = 0; k < k_fold; k++ )
1847 {
1848 int b0 = ((k+1)*n0 + k_fold/2)/k_fold, b1 = ((k+1)*n1 + k_fold/2)/k_fold;
1849 int a = (int)sidx.size(), b = a + (b0 - a0) + (b1 - a1);
1850 for( int i = a0; i < b0; i++ )
1851 sidx.push_back(sidx0[i]);
1852 for( int i = a1; i < b1; i++ )
1853 sidx.push_back(sidx1[i]);
1854 for( int i = 0; i < (b - a); i++ )
1855 {
1856 int i1 = rng.uniform(a, b);
1857 int i2 = rng.uniform(a, b);
1858 std::swap(sidx[i1], sidx[i2]);
1859 }
1860 a0 = b0; a1 = b1;
1861 }
1862 }
1863
1864 // If grid.minVal == grid.maxVal, this will allow one and only one pass through the loop with params.var = grid.minVal.
1865 #define FOR_IN_GRID(var, grid) \
1866 for( params.var = grid.minVal; params.var == grid.minVal || params.var < grid.maxVal; params.var = (grid.minVal == grid.maxVal) ? grid.maxVal + 1 : params.var * grid.logStep )
1867
1868 // Create the list of parameters to test
1869 std::vector<SvmParams> parameters;
1870 FOR_IN_GRID(C, C_grid)
1871 FOR_IN_GRID(gamma, gamma_grid)
1872 FOR_IN_GRID(p, p_grid)
1873 FOR_IN_GRID(nu, nu_grid)
1874 FOR_IN_GRID(coef0, coef_grid)
1875 FOR_IN_GRID(degree, degree_grid)
1876 {
1877 parameters.push_back(params);
1878 }
1879
1880 std::vector<double> result(parameters.size());
1881 TrainAutoBody invoker(parameters, samples, responses, class_labels, sidx,
1882 is_classification, k_fold, result);
1883 parallel_for_(cv::Range(0,(int)parameters.size()), invoker);
1884
1885 // Extract the best parameters
1886 SvmParams best_params = params;
1887 double min_error = FLT_MAX;
1888 for( int i = 0; i < (int)result.size(); i++ )
1889 {
1890 if( result[i] < min_error )
1891 {
1892 min_error = result[i];
1893 best_params = parameters[i];
1894 }
1895 }
1896
1897 class_labels = class_labels0;
1898 setParams(best_params);
1899 return do_train( samples, responses );
1900 }
1901
1902 struct PredictBody : ParallelLoopBody
1903 {
PredictBodycv::ml::CV_FINAL::PredictBody1904 PredictBody( const SVMImpl* _svm, const Mat& _samples, Mat& _results, bool _returnDFVal )
1905 {
1906 svm = _svm;
1907 results = &_results;
1908 samples = &_samples;
1909 returnDFVal = _returnDFVal;
1910 }
1911
operator ()cv::ml::CV_FINAL::PredictBody1912 void operator()(const Range& range) const CV_OVERRIDE
1913 {
1914 int svmType = svm->params.svmType;
1915 int sv_total = svm->sv.rows;
1916 int class_count = !svm->class_labels.empty() ? (int)svm->class_labels.total() : svmType == ONE_CLASS ? 1 : 0;
1917
1918 AutoBuffer<float> _buffer(sv_total + (class_count+1)*2);
1919 float* buffer = _buffer.data();
1920
1921 int i, j, dfi, k, si;
1922
1923 if( svmType == EPS_SVR || svmType == NU_SVR || svmType == ONE_CLASS )
1924 {
1925 for( si = range.start; si < range.end; si++ )
1926 {
1927 const float* row_sample = samples->ptr<float>(si);
1928 svm->kernel->calc( sv_total, svm->var_count, svm->sv.ptr<float>(), row_sample, buffer );
1929
1930 const SVMImpl::DecisionFunc* df = &svm->decision_func[0];
1931 double sum = -df->rho;
1932 for( i = 0; i < sv_total; i++ )
1933 sum += buffer[i]*svm->df_alpha[i];
1934 float result = svm->params.svmType == ONE_CLASS && !returnDFVal ? (float)(sum > 0) : (float)sum;
1935 results->at<float>(si) = result;
1936 }
1937 }
1938 else if( svmType == C_SVC || svmType == NU_SVC )
1939 {
1940 int* vote = (int*)(buffer + sv_total);
1941
1942 for( si = range.start; si < range.end; si++ )
1943 {
1944 svm->kernel->calc( sv_total, svm->var_count, svm->sv.ptr<float>(),
1945 samples->ptr<float>(si), buffer );
1946 double sum = 0.;
1947
1948 memset( vote, 0, class_count*sizeof(vote[0]));
1949
1950 for( i = dfi = 0; i < class_count; i++ )
1951 {
1952 for( j = i+1; j < class_count; j++, dfi++ )
1953 {
1954 const DecisionFunc& df = svm->decision_func[dfi];
1955 sum = -df.rho;
1956 int sv_count = svm->getSVCount(dfi);
1957 CV_DbgAssert(sv_count > 0);
1958 const double* alpha = &svm->df_alpha[df.ofs];
1959 const int* sv_index = &svm->df_index[df.ofs];
1960 for( k = 0; k < sv_count; k++ )
1961 sum += alpha[k]*buffer[sv_index[k]];
1962
1963 vote[sum > 0 ? i : j]++;
1964 }
1965 }
1966
1967 for( i = 1, k = 0; i < class_count; i++ )
1968 {
1969 if( vote[i] > vote[k] )
1970 k = i;
1971 }
1972 float result = returnDFVal && class_count == 2 ?
1973 (float)sum : (float)(svm->class_labels.at<int>(k));
1974 results->at<float>(si) = result;
1975 }
1976 }
1977 else
1978 CV_Error( CV_StsBadArg, "INTERNAL ERROR: Unknown SVM type, "
1979 "the SVM structure is probably corrupted" );
1980 }
1981
1982 const SVMImpl* svm;
1983 const Mat* samples;
1984 Mat* results;
1985 bool returnDFVal;
1986 };
1987
trainAuto(InputArray samples,int layout,InputArray responses,int kfold,Ptr<ParamGrid> Cgrid,Ptr<ParamGrid> gammaGrid,Ptr<ParamGrid> pGrid,Ptr<ParamGrid> nuGrid,Ptr<ParamGrid> coeffGrid,Ptr<ParamGrid> degreeGrid,bool balanced)1988 bool trainAuto(InputArray samples, int layout,
1989 InputArray responses, int kfold, Ptr<ParamGrid> Cgrid,
1990 Ptr<ParamGrid> gammaGrid, Ptr<ParamGrid> pGrid, Ptr<ParamGrid> nuGrid,
1991 Ptr<ParamGrid> coeffGrid, Ptr<ParamGrid> degreeGrid, bool balanced) CV_OVERRIDE
1992 {
1993 Ptr<TrainData> data = TrainData::create(samples, layout, responses);
1994 return this->trainAuto(
1995 data, kfold,
1996 *Cgrid.get(),
1997 *gammaGrid.get(),
1998 *pGrid.get(),
1999 *nuGrid.get(),
2000 *coeffGrid.get(),
2001 *degreeGrid.get(),
2002 balanced);
2003 }
2004
2005
predict(InputArray _samples,OutputArray _results,int flags) const2006 float predict( InputArray _samples, OutputArray _results, int flags ) const CV_OVERRIDE
2007 {
2008 float result = 0;
2009 Mat samples = _samples.getMat(), results;
2010 int nsamples = samples.rows;
2011 bool returnDFVal = (flags & RAW_OUTPUT) != 0;
2012
2013 CV_Assert( samples.cols == var_count && samples.type() == CV_32F );
2014
2015 if( _results.needed() )
2016 {
2017 _results.create( nsamples, 1, samples.type() );
2018 results = _results.getMat();
2019 }
2020 else
2021 {
2022 CV_Assert( nsamples == 1 );
2023 results = Mat(1, 1, CV_32F, &result);
2024 }
2025
2026 PredictBody invoker(this, samples, results, returnDFVal);
2027 if( nsamples < 10 )
2028 invoker(Range(0, nsamples));
2029 else
2030 parallel_for_(Range(0, nsamples), invoker);
2031 return result;
2032 }
2033
getDecisionFunction(int i,OutputArray _alpha,OutputArray _svidx) const2034 double getDecisionFunction(int i, OutputArray _alpha, OutputArray _svidx ) const CV_OVERRIDE
2035 {
2036 CV_Assert( 0 <= i && i < (int)decision_func.size());
2037 const DecisionFunc& df = decision_func[i];
2038 int count = getSVCount(i);
2039 Mat(1, count, CV_64F, (double*)&df_alpha[df.ofs]).copyTo(_alpha);
2040 Mat(1, count, CV_32S, (int*)&df_index[df.ofs]).copyTo(_svidx);
2041 return df.rho;
2042 }
2043
write_params(FileStorage & fs) const2044 void write_params( FileStorage& fs ) const
2045 {
2046 int svmType = params.svmType;
2047 int kernelType = params.kernelType;
2048
2049 String svm_type_str =
2050 svmType == C_SVC ? "C_SVC" :
2051 svmType == NU_SVC ? "NU_SVC" :
2052 svmType == ONE_CLASS ? "ONE_CLASS" :
2053 svmType == EPS_SVR ? "EPS_SVR" :
2054 svmType == NU_SVR ? "NU_SVR" : format("Unknown_%d", svmType);
2055 String kernel_type_str =
2056 kernelType == LINEAR ? "LINEAR" :
2057 kernelType == POLY ? "POLY" :
2058 kernelType == RBF ? "RBF" :
2059 kernelType == SIGMOID ? "SIGMOID" :
2060 kernelType == CHI2 ? "CHI2" :
2061 kernelType == INTER ? "INTER" : format("Unknown_%d", kernelType);
2062
2063 fs << "svmType" << svm_type_str;
2064
2065 // save kernel
2066 fs << "kernel" << "{" << "type" << kernel_type_str;
2067
2068 if( kernelType == POLY )
2069 fs << "degree" << params.degree;
2070
2071 if( kernelType != LINEAR )
2072 fs << "gamma" << params.gamma;
2073
2074 if( kernelType == POLY || kernelType == SIGMOID )
2075 fs << "coef0" << params.coef0;
2076
2077 fs << "}";
2078
2079 if( svmType == C_SVC || svmType == EPS_SVR || svmType == NU_SVR )
2080 fs << "C" << params.C;
2081
2082 if( svmType == NU_SVC || svmType == ONE_CLASS || svmType == NU_SVR )
2083 fs << "nu" << params.nu;
2084
2085 if( svmType == EPS_SVR )
2086 fs << "p" << params.p;
2087
2088 fs << "term_criteria" << "{:";
2089 if( params.termCrit.type & TermCriteria::EPS )
2090 fs << "epsilon" << params.termCrit.epsilon;
2091 if( params.termCrit.type & TermCriteria::COUNT )
2092 fs << "iterations" << params.termCrit.maxCount;
2093 fs << "}";
2094 }
2095
isTrained() const2096 bool isTrained() const CV_OVERRIDE
2097 {
2098 return !sv.empty();
2099 }
2100
isClassifier() const2101 bool isClassifier() const CV_OVERRIDE
2102 {
2103 return params.svmType == C_SVC || params.svmType == NU_SVC || params.svmType == ONE_CLASS;
2104 }
2105
getVarCount() const2106 int getVarCount() const CV_OVERRIDE
2107 {
2108 return var_count;
2109 }
2110
getDefaultName() const2111 String getDefaultName() const CV_OVERRIDE
2112 {
2113 return "opencv_ml_svm";
2114 }
2115
write(FileStorage & fs) const2116 void write( FileStorage& fs ) const CV_OVERRIDE
2117 {
2118 int class_count = !class_labels.empty() ? (int)class_labels.total() :
2119 params.svmType == ONE_CLASS ? 1 : 0;
2120 if( !isTrained() )
2121 CV_Error( CV_StsParseError, "SVM model data is invalid, check sv_count, var_* and class_count tags" );
2122
2123 writeFormat(fs);
2124 write_params( fs );
2125
2126 fs << "var_count" << var_count;
2127
2128 if( class_count > 0 )
2129 {
2130 fs << "class_count" << class_count;
2131
2132 if( !class_labels.empty() )
2133 fs << "class_labels" << class_labels;
2134
2135 if( !params.classWeights.empty() )
2136 fs << "class_weights" << params.classWeights;
2137 }
2138
2139 // write the joint collection of support vectors
2140 int i, sv_total = sv.rows;
2141 fs << "sv_total" << sv_total;
2142 fs << "support_vectors" << "[";
2143 for( i = 0; i < sv_total; i++ )
2144 {
2145 fs << "[:";
2146 fs.writeRaw("f", sv.ptr(i), sv.cols*sv.elemSize());
2147 fs << "]";
2148 }
2149 fs << "]";
2150
2151 if ( !uncompressed_sv.empty() )
2152 {
2153 // write the joint collection of uncompressed support vectors
2154 int uncompressed_sv_total = uncompressed_sv.rows;
2155 fs << "uncompressed_sv_total" << uncompressed_sv_total;
2156 fs << "uncompressed_support_vectors" << "[";
2157 for( i = 0; i < uncompressed_sv_total; i++ )
2158 {
2159 fs << "[:";
2160 fs.writeRaw("f", uncompressed_sv.ptr(i), uncompressed_sv.cols*uncompressed_sv.elemSize());
2161 fs << "]";
2162 }
2163 fs << "]";
2164 }
2165
2166 // write decision functions
2167 int df_count = (int)decision_func.size();
2168
2169 fs << "decision_functions" << "[";
2170 for( i = 0; i < df_count; i++ )
2171 {
2172 const DecisionFunc& df = decision_func[i];
2173 int sv_count = getSVCount(i);
2174 fs << "{" << "sv_count" << sv_count
2175 << "rho" << df.rho
2176 << "alpha" << "[:";
2177 fs.writeRaw("d", (const uchar*)&df_alpha[df.ofs], sv_count*sizeof(df_alpha[0]));
2178 fs << "]";
2179 if( class_count >= 2 )
2180 {
2181 fs << "index" << "[:";
2182 fs.writeRaw("i", (const uchar*)&df_index[df.ofs], sv_count*sizeof(df_index[0]));
2183 fs << "]";
2184 }
2185 else
2186 CV_Assert( sv_count == sv_total );
2187 fs << "}";
2188 }
2189 fs << "]";
2190 }
2191
read_params(const FileNode & fn)2192 void read_params( const FileNode& fn )
2193 {
2194 SvmParams _params;
2195
2196 // check for old naming
2197 String svm_type_str = (String)(fn["svm_type"].empty() ? fn["svmType"] : fn["svm_type"]);
2198 int svmType =
2199 svm_type_str == "C_SVC" ? C_SVC :
2200 svm_type_str == "NU_SVC" ? NU_SVC :
2201 svm_type_str == "ONE_CLASS" ? ONE_CLASS :
2202 svm_type_str == "EPS_SVR" ? EPS_SVR :
2203 svm_type_str == "NU_SVR" ? NU_SVR : -1;
2204
2205 if( svmType < 0 )
2206 CV_Error( CV_StsParseError, "Missing or invalid SVM type" );
2207
2208 FileNode kernel_node = fn["kernel"];
2209 if( kernel_node.empty() )
2210 CV_Error( CV_StsParseError, "SVM kernel tag is not found" );
2211
2212 String kernel_type_str = (String)kernel_node["type"];
2213 int kernelType =
2214 kernel_type_str == "LINEAR" ? LINEAR :
2215 kernel_type_str == "POLY" ? POLY :
2216 kernel_type_str == "RBF" ? RBF :
2217 kernel_type_str == "SIGMOID" ? SIGMOID :
2218 kernel_type_str == "CHI2" ? CHI2 :
2219 kernel_type_str == "INTER" ? INTER : CUSTOM;
2220
2221 if( kernelType == CUSTOM )
2222 CV_Error( CV_StsParseError, "Invalid SVM kernel type (or custom kernel)" );
2223
2224 _params.svmType = svmType;
2225 _params.kernelType = kernelType;
2226 _params.degree = (double)kernel_node["degree"];
2227 _params.gamma = (double)kernel_node["gamma"];
2228 _params.coef0 = (double)kernel_node["coef0"];
2229
2230 _params.C = (double)fn["C"];
2231 _params.nu = (double)fn["nu"];
2232 _params.p = (double)fn["p"];
2233 _params.classWeights = Mat();
2234
2235 FileNode tcnode = fn["term_criteria"];
2236 if( !tcnode.empty() )
2237 {
2238 _params.termCrit.epsilon = (double)tcnode["epsilon"];
2239 _params.termCrit.maxCount = (int)tcnode["iterations"];
2240 _params.termCrit.type = (_params.termCrit.epsilon > 0 ? TermCriteria::EPS : 0) +
2241 (_params.termCrit.maxCount > 0 ? TermCriteria::COUNT : 0);
2242 }
2243 else
2244 _params.termCrit = TermCriteria( TermCriteria::EPS + TermCriteria::COUNT, 1000, FLT_EPSILON );
2245
2246 setParams( _params );
2247 }
2248
read(const FileNode & fn)2249 void read( const FileNode& fn ) CV_OVERRIDE
2250 {
2251 clear();
2252
2253 // read SVM parameters
2254 read_params( fn );
2255
2256 // and top-level data
2257 int i, sv_total = (int)fn["sv_total"];
2258 var_count = (int)fn["var_count"];
2259 int class_count = (int)fn["class_count"];
2260
2261 if( sv_total <= 0 || var_count <= 0 )
2262 CV_Error( CV_StsParseError, "SVM model data is invalid, check sv_count, var_* and class_count tags" );
2263
2264 FileNode m = fn["class_labels"];
2265 if( !m.empty() )
2266 m >> class_labels;
2267 m = fn["class_weights"];
2268 if( !m.empty() )
2269 m >> params.classWeights;
2270
2271 if( class_count > 1 && (class_labels.empty() || (int)class_labels.total() != class_count))
2272 CV_Error( CV_StsParseError, "Array of class labels is missing or invalid" );
2273
2274 // read support vectors
2275 FileNode sv_node = fn["support_vectors"];
2276
2277 CV_Assert((int)sv_node.size() == sv_total);
2278
2279 sv.create(sv_total, var_count, CV_32F);
2280 FileNodeIterator sv_it = sv_node.begin();
2281 for( i = 0; i < sv_total; i++, ++sv_it )
2282 {
2283 (*sv_it).readRaw("f", sv.ptr(i), var_count*sv.elemSize());
2284 }
2285
2286 int uncompressed_sv_total = (int)fn["uncompressed_sv_total"];
2287
2288 if( uncompressed_sv_total > 0 )
2289 {
2290 // read uncompressed support vectors
2291 FileNode uncompressed_sv_node = fn["uncompressed_support_vectors"];
2292
2293 CV_Assert((int)uncompressed_sv_node.size() == uncompressed_sv_total);
2294 uncompressed_sv.create(uncompressed_sv_total, var_count, CV_32F);
2295
2296 FileNodeIterator uncompressed_sv_it = uncompressed_sv_node.begin();
2297 for( i = 0; i < uncompressed_sv_total; i++, ++uncompressed_sv_it )
2298 {
2299 (*uncompressed_sv_it).readRaw("f", uncompressed_sv.ptr(i), var_count*uncompressed_sv.elemSize());
2300 }
2301 }
2302
2303 // read decision functions
2304 int df_count = class_count > 1 ? class_count*(class_count-1)/2 : 1;
2305 FileNode df_node = fn["decision_functions"];
2306
2307 CV_Assert((int)df_node.size() == df_count);
2308
2309 FileNodeIterator df_it = df_node.begin();
2310 for( i = 0; i < df_count; i++, ++df_it )
2311 {
2312 FileNode dfi = *df_it;
2313 DecisionFunc df;
2314 int sv_count = (int)dfi["sv_count"];
2315 int ofs = (int)df_index.size();
2316 df.rho = (double)dfi["rho"];
2317 df.ofs = ofs;
2318 df_index.resize(ofs + sv_count);
2319 df_alpha.resize(ofs + sv_count);
2320 dfi["alpha"].readRaw("d", (uchar*)&df_alpha[ofs], sv_count*sizeof(df_alpha[0]));
2321 if( class_count >= 2 )
2322 dfi["index"].readRaw("i", (uchar*)&df_index[ofs], sv_count*sizeof(df_index[0]));
2323 decision_func.push_back(df);
2324 }
2325 if( class_count < 2 )
2326 setRangeVector(df_index, sv_total);
2327 if( (int)fn["optimize_linear"] != 0 )
2328 optimize_linear_svm();
2329 }
2330
2331 SvmParams params;
2332 Mat class_labels;
2333 int var_count;
2334 Mat sv, uncompressed_sv;
2335 vector<DecisionFunc> decision_func;
2336 vector<double> df_alpha;
2337 vector<int> df_index;
2338
2339 Ptr<Kernel> kernel;
2340 };
2341
2342
create()2343 Ptr<SVM> SVM::create()
2344 {
2345 return makePtr<SVMImpl>();
2346 }
2347
load(const String & filepath)2348 Ptr<SVM> SVM::load(const String& filepath)
2349 {
2350 FileStorage fs;
2351 fs.open(filepath, FileStorage::READ);
2352
2353 Ptr<SVM> svm = makePtr<SVMImpl>();
2354
2355 ((SVMImpl*)svm.get())->read(fs.getFirstTopLevelNode());
2356 return svm;
2357 }
2358
2359
2360 }
2361 }
2362
2363 /* End of file. */
2364