{

if( pg.minVal > pg.maxVal )

CV_Error( CV_StsBadArg, "Lower bound of the grid must be less then the upper one" );

if( pg.minVal < DBL_EPSILON )

CV_Error( CV_StsBadArg, "Lower bound of the grid must be positive" );

if( pg.logStep < 1. + FLT_EPSILON )

CV_Error( CV_StsBadArg, "Grid step must greater then 1" );

{

int svmType;

int kernelType;

double gamma;

double coef0;

double degree;

double C;

double nu;

double p;

Mat classWeights;

TermCriteria termCrit;

SvmParams()

{

svmType = SVM::C_SVC;

kernelType = SVM::RBF;

degree = 0;

gamma = 1;

coef0 = 0;

C = 1;

nu = 0;

p = 0;

termCrit = TermCriteria( CV_TERMCRIT_ITER+CV_TERMCRIT_EPS, 1000, FLT_EPSILON );

}

SvmParams( int _svmType, int _kernelType,

double _degree, double _gamma, double _coef0,

double _Con, double _nu, double _p,

const Mat& _classWeights, TermCriteria _termCrit )

{

svmType = _svmType;

kernelType = _kernelType;

degree = _degree;

gamma = _gamma;

coef0 = _coef0;

C = _Con;

nu = _nu;

p = _p;

classWeights = _classWeights;

termCrit = _termCrit;

}

{

SVMKernelImpl( const SvmParams& _params = SvmParams() )

{

params = _params;

}

int getType() const

{

return params.kernelType;

}

void calc_non_rbf_base( int vcount, int var_count, const float* vecs,

const float* another, Qfloat* results,

double alpha, double beta )

{

int j, k;

for( j = 0; j < vcount; j++ )

{

const float* sample = &vecs[j*var_count];

double s = 0;

for( k = 0; k <= var_count - 4; k += 4 )

s += sample[k]*another[k] + sample[k+1]*another[k+1] +

sample[k+2]*another[k+2] + sample[k+3]*another[k+3];

for( ; k < var_count; k++ )

s += sample[k]*another[k];

results[j] = (Qfloat)(s*alpha + beta);

}

}

void calc_linear( int vcount, int var_count, const float* vecs,

const float* another, Qfloat* results )

{

calc_non_rbf_base( vcount, var_count, vecs, another, results, 1, 0 );

}

void calc_poly( int vcount, int var_count, const float* vecs,

const float* another, Qfloat* results )

{

Mat R( 1, vcount, QFLOAT_TYPE, results );

calc_non_rbf_base( vcount, var_count, vecs, another, results, params.gamma, params.coef0 );

if( vcount > 0 )

pow( R, params.degree, R );

}

void calc_sigmoid( int vcount, int var_count, const float* vecs,

const float* another, Qfloat* results )

{

int j;

calc_non_rbf_base( vcount, var_count, vecs, another, results,

-2*params.gamma, -2*params.coef0 );

// TODO: speedup this

for( j = 0; j < vcount; j++ )

{

Qfloat t = results[j];

Qfloat e = std::exp(-std::abs(t));

if( t > 0 )

results[j] = (Qfloat)((1. - e)/(1. + e));

else

results[j] = (Qfloat)((e - 1.)/(e + 1.));

}

}

void calc_rbf( int vcount, int var_count, const float* vecs,

const float* another, Qfloat* results )

{

double gamma = -params.gamma;

int j, k;

for( j = 0; j < vcount; j++ )

{

const float* sample = &vecs[j*var_count];

double s = 0;

for( k = 0; k <= var_count - 4; k += 4 )

{

double t0 = sample[k] - another[k];

double t1 = sample[k+1] - another[k+1];

s += t0*t0 + t1*t1;

t0 = sample[k+2] - another[k+2];

t1 = sample[k+3] - another[k+3];

s += t0*t0 + t1*t1;

}

for( ; k < var_count; k++ )

{

double t0 = sample[k] - another[k];

s += t0*t0;

}

results[j] = (Qfloat)(s*gamma);

}

if( vcount > 0 )

{

Mat R( 1, vcount, QFLOAT_TYPE, results );

exp( R, R );

}

}

/// Histogram intersection kernel

void calc_intersec( int vcount, int var_count, const float* vecs,

const float* another, Qfloat* results )

{

int j, k;

for( j = 0; j < vcount; j++ )

{

const float* sample = &vecs[j*var_count];

double s = 0;

for( k = 0; k <= var_count - 4; k += 4 )

s += std::min(sample[k],another[k]) + std::min(sample[k+1],another[k+1]) +

std::min(sample[k+2],another[k+2]) + std::min(sample[k+3],another[k+3]);

for( ; k < var_count; k++ )

s += std::min(sample[k],another[k]);

results[j] = (Qfloat)(s);

}

}

/// Exponential chi2 kernel

void calc_chi2( int vcount, int var_count, const float* vecs,

const float* another, Qfloat* results )

{

Mat R( 1, vcount, QFLOAT_TYPE, results );

double gamma = -params.gamma;

int j, k;

for( j = 0; j < vcount; j++ )

{

const float* sample = &vecs[j*var_count];

double chi2 = 0;

for(k = 0 ; k < var_count; k++ )

{

double d = sample[k]-another[k];

double devisor = sample[k]+another[k];

/// if devisor == 0, the Chi2 distance would be zero,

// but calculation would rise an error because of deviding by zero

if (devisor != 0)

{

chi2 += d*d/devisor;

}

}

results[j] = (Qfloat) (gamma*chi2);

}

if( vcount > 0 )

exp( R, R );

}

void calc( int vcount, int var_count, const float* vecs,

const float* another, Qfloat* results )

{

switch( params.kernelType )

{

case SVM::LINEAR:

calc_linear(vcount, var_count, vecs, another, results);

break;

case SVM::RBF:

calc_rbf(vcount, var_count, vecs, another, results);

break;

case SVM::POLY:

calc_poly(vcount, var_count, vecs, another, results);

break;

case SVM::SIGMOID:

calc_sigmoid(vcount, var_count, vecs, another, results);

break;

case SVM::CHI2:

calc_chi2(vcount, var_count, vecs, another, results);

break;

case SVM::INTER:

calc_intersec(vcount, var_count, vecs, another, results);

break;

default:

CV_Error(CV_StsBadArg, "Unknown kernel type");

}

const Qfloat max_val = (Qfloat)(FLT_MAX*1e-3);

for( int j = 0; j < vcount; j++ )

{

if( results[j] > max_val )

results[j] = max_val;

}

}

SvmParams params;

vector<int>& sidx_all, vector<int>& class_ranges )

{

int i, nsamples = _samples.rows;

CV_Assert( _responses.isContinuous() && _responses.checkVector(1, CV_32S) == nsamples );

setRangeVector(sidx_all, nsamples);

const int* rptr = _responses.ptr<int>();

std::sort(sidx_all.begin(), sidx_all.end(), cmp_lt_idx<int>(rptr));

class_ranges.clear();

class_ranges.push_back(0);

for( i = 0; i < nsamples; i++ )

{

if( i == nsamples-1 || rptr[sidx_all[i]] != rptr[sidx_all[i+1]] )

class_ranges.push_back(i+1);

}

{

ParamGrid grid = getDefaultGrid(param_id); // this is not a nice solution..

return makePtr<ParamGrid>(grid.minVal, grid.maxVal, grid.logStep);

{

ParamGrid grid;

if( param_id == SVM::C )

{

grid.minVal = 0.1;

grid.maxVal = 500;

grid.logStep = 5; // total iterations = 5

}

else if( param_id == SVM::GAMMA )

{

grid.minVal = 1e-5;

grid.maxVal = 0.6;

grid.logStep = 15; // total iterations = 4

}

else if( param_id == SVM::P )

{

grid.minVal = 0.01;

grid.maxVal = 100;

grid.logStep = 7; // total iterations = 4

}

else if( param_id == SVM::NU )

{

grid.minVal = 0.01;

grid.maxVal = 0.2;

grid.logStep = 3; // total iterations = 3

}

else if( param_id == SVM::COEF )

{

grid.minVal = 0.1;

grid.maxVal = 300;

grid.logStep = 14; // total iterations = 3

}

else if( param_id == SVM::DEGREE )

{

grid.minVal = 0.01;

grid.maxVal = 4;

grid.logStep = 7; // total iterations = 3

}

else

cvError( CV_StsBadArg, "SVM::getDefaultGrid", "Invalid type of parameter "

"(use one of SVM::C, SVM::GAMMA et al.)", __FILE__, __LINE__ );

return grid;

{

struct DecisionFunc

{

DecisionFunc(double _rho, int _ofs) : rho(_rho), ofs(_ofs) {}

DecisionFunc() : rho(0.), ofs(0) {}

double rho;

int ofs;

};

// Generalized SMO+SVMlight algorithm

// Solves:

//

// min [0.5(\alpha^T Q \alpha) + b^T \alpha]

//

// y^T \alpha = \delta

// y_i = +1 or -1

// 0 <= alpha_i <= Cp for y_i = 1

// 0 <= alpha_i <= Cn for y_i = -1

//

// Given:

//

// Q, b, y, Cp, Cn, and an initial feasible point \alpha

// l is the size of vectors and matrices

// eps is the stopping criterion

//

// solution will be put in \alpha, objective value will be put in obj

//

class Solver

{

public:

enum { MIN_CACHE_SIZE = (40 << 20) /* 40Mb */, MAX_CACHE_SIZE = (500 << 20) /* 500Mb */ };

typedef bool (Solver::*SelectWorkingSet)( int& i, int& j );

typedef Qfloat* (Solver::*GetRow)( int i, Qfloat* row, Qfloat* dst, bool existed );

typedef void (Solver::*CalcRho)( double& rho, double& r );

struct KernelRow

{

KernelRow() { idx = -1; prev = next = 0; }

KernelRow(int _idx, int _prev, int _next) : idx(_idx), prev(_prev), next(_next) {}

int idx;

int prev;

int next;

};

struct SolutionInfo

{

SolutionInfo() { obj = rho = upper_bound_p = upper_bound_n = r = 0; }

double obj;

double rho;

double upper_bound_p;

double upper_bound_n;

double r; // for Solver_NU

};

void clear()

{

alpha_vec = 0;

select_working_set_func = 0;

calc_rho_func = 0;

get_row_func = 0;

lru_cache.clear();

}

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 )

{

clear();

samples = _samples;

sample_count = samples.rows;

var_count = samples.cols;

y_vec = _y;

alpha_vec = &_alpha;

alpha_count = (int)alpha_vec->size();

b_vec = _b;

kernel = _kernel;

C[0] = _Cn;

C[1] = _Cp;

eps = _termCrit.epsilon;

max_iter = _termCrit.maxCount;

G_vec.resize(alpha_count);

alpha_status_vec.resize(alpha_count);

buf[0].resize(sample_count*2);

buf[1].resize(sample_count*2);

select_working_set_func = _select_working_set;

CV_Assert(select_working_set_func != 0);

calc_rho_func = _calc_rho;

CV_Assert(calc_rho_func != 0);

get_row_func = _get_row;

CV_Assert(get_row_func != 0);

// assume that for large training sets ~25% of Q matrix is used

int64 csize = (int64)sample_count*sample_count/4;

csize = std::max(csize, (int64)(MIN_CACHE_SIZE/sizeof(Qfloat)) );

csize = std::min(csize, (int64)(MAX_CACHE_SIZE/sizeof(Qfloat)) );

max_cache_size = (int)((csize + sample_count-1)/sample_count);

max_cache_size = std::min(std::max(max_cache_size, 1), sample_count);

cache_size = 0;

lru_cache.clear();

lru_cache.resize(sample_count+1, KernelRow(-1, 0, 0));

lru_first = lru_last = 0;

lru_cache_data.create(max_cache_size, sample_count, QFLOAT_TYPE);

}

Qfloat* get_row_base( int i, bool* _existed )

{

int i1 = i < sample_count ? i : i - sample_count;

KernelRow& kr = lru_cache[i1+1];

if( _existed )

*_existed = kr.idx >= 0;

if( kr.idx < 0 )

{

if( cache_size < max_cache_size )

{

kr.idx = cache_size;

cache_size++;

if (!lru_last)

lru_last = i1+1;

}

else

{

KernelRow& last = lru_cache[lru_last];

kr.idx = last.idx;

last.idx = -1;

lru_cache[last.prev].next = 0;

lru_last = last.prev;

last.prev = 0;

last.next = 0;

}

kernel->calc( sample_count, var_count, samples.ptr<float>(),

samples.ptr<float>(i1), lru_cache_data.ptr<Qfloat>(kr.idx) );

}

else

{

if( kr.next )

lru_cache[kr.next].prev = kr.prev;

else

lru_last = kr.prev;

if( kr.prev )

lru_cache[kr.prev].next = kr.next;

else

lru_first = kr.next;

}

if (lru_first)

lru_cache[lru_first].prev = i1+1;

kr.next = lru_first;

kr.prev = 0;

lru_first = i1+1;

return lru_cache_data.ptr<Qfloat>(kr.idx);

}

Qfloat* get_row_svc( int i, Qfloat* row, Qfloat*, bool existed )

{

if( !existed )

{

const schar* _y = &y_vec[0];

int j, len = sample_count;

if( _y[i] > 0 )

{

for( j = 0; j < len; j++ )

row[j] = _y[j]*row[j];

}

else

{

for( j = 0; j < len; j++ )

row[j] = -_y[j]*row[j];

}

}

return row;

}

Qfloat* get_row_one_class( int, Qfloat* row, Qfloat*, bool )

{

return row;

}

Qfloat* get_row_svr( int i, Qfloat* row, Qfloat* dst, bool )

{

int j, len = sample_count;

Qfloat* dst_pos = dst;

Qfloat* dst_neg = dst + len;

if( i >= len )

std::swap(dst_pos, dst_neg);

for( j = 0; j < len; j++ )

{

Qfloat t = row[j];

dst_pos[j] = t;

dst_neg[j] = -t;

}

return dst;

}

Qfloat* get_row( int i, float* dst )

{

bool existed = false;

float* row = get_row_base( i, &existed );

return (this->*get_row_func)( i, row, dst, existed );

}

#undef is_upper_bound

#define is_upper_bound(i) (alpha_status[i] > 0)

#undef is_lower_bound

#define is_lower_bound(i) (alpha_status[i] < 0)

#undef is_free

#define is_free(i) (alpha_status[i] == 0)

#undef get_C

#define get_C(i) (C[y[i]>0])

#undef update_alpha_status

#define update_alpha_status(i) \

alpha_status[i] = (schar)(alpha[i] >= get_C(i) ? 1 : alpha[i] <= 0 ? -1 : 0)

#undef reconstruct_gradient

#define reconstruct_gradient() /* empty for now */

bool solve_generic( SolutionInfo& si )

{

const schar* y = &y_vec[0];

double* alpha = &alpha_vec->at(0);

schar* alpha_status = &alpha_status_vec[0];

double* G = &G_vec[0];

double* b = &b_vec[0];

int iter = 0;

int i, j, k;

// 1. initialize gradient and alpha status

for( i = 0; i < alpha_count; i++ )

{

update_alpha_status(i);

G[i] = b[i];

if( fabs(G[i]) > 1e200 )

return false;

}

for( i = 0; i < alpha_count; i++ )

{

if( !is_lower_bound(i) )

{

const Qfloat *Q_i = get_row( i, &buf[0][0] );

double alpha_i = alpha[i];

for( j = 0; j < alpha_count; j++ )

G[j] += alpha_i*Q_i[j];

}

}

// 2. optimization loop

for(;;)

{

const Qfloat *Q_i, *Q_j;

double C_i, C_j;

double old_alpha_i, old_alpha_j, alpha_i, alpha_j;

double delta_alpha_i, delta_alpha_j;

#ifdef _DEBUG

for( i = 0; i < alpha_count; i++ )

{

if( fabs(G[i]) > 1e+300 )

return false;

if( fabs(alpha[i]) > 1e16 )

return false;

}

#endif

if( (this->*select_working_set_func)( i, j ) != 0 || iter++ >= max_iter )

break;

Q_i = get_row( i, &buf[0][0] );

Q_j = get_row( j, &buf[1][0] );

C_i = get_C(i);

C_j = get_C(j);

alpha_i = old_alpha_i = alpha[i];

alpha_j = old_alpha_j = alpha[j];

if( y[i] != y[j] )

{

double denom = Q_i[i]+Q_j[j]+2*Q_i[j];

double delta = (-G[i]-G[j])/MAX(fabs(denom),FLT_EPSILON);

double diff = alpha_i - alpha_j;

alpha_i += delta;

alpha_j += delta;

if( diff > 0 && alpha_j < 0 )

{

alpha_j = 0;

alpha_i = diff;

}

else if( diff <= 0 && alpha_i < 0 )

{

alpha_i = 0;

alpha_j = -diff;

}

if( diff > C_i - C_j && alpha_i > C_i )

{

alpha_i = C_i;

alpha_j = C_i - diff;

}

else if( diff <= C_i - C_j && alpha_j > C_j )

{

alpha_j = C_j;

alpha_i = C_j + diff;

}

}

else

{

double denom = Q_i[i]+Q_j[j]-2*Q_i[j];

double delta = (G[i]-G[j])/MAX(fabs(denom),FLT_EPSILON);

double sum = alpha_i + alpha_j;

alpha_i -= delta;

alpha_j += delta;

if( sum > C_i && alpha_i > C_i )

{

alpha_i = C_i;

alpha_j = sum - C_i;

}

else if( sum <= C_i && alpha_j < 0)

{

alpha_j = 0;

alpha_i = sum;

}

if( sum > C_j && alpha_j > C_j )

{

alpha_j = C_j;

alpha_i = sum - C_j;

}

else if( sum <= C_j && alpha_i < 0 )

{

alpha_i = 0;

alpha_j = sum;

}

}

// update alpha

alpha[i] = alpha_i;

alpha[j] = alpha_j;

update_alpha_status(i);

update_alpha_status(j);

// update G

delta_alpha_i = alpha_i - old_alpha_i;

delta_alpha_j = alpha_j - old_alpha_j;

for( k = 0; k < alpha_count; k++ )

G[k] += Q_i[k]*delta_alpha_i + Q_j[k]*delta_alpha_j;

}

// calculate rho

(this->*calc_rho_func)( si.rho, si.r );

// calculate objective value

for( i = 0, si.obj = 0; i < alpha_count; i++ )

si.obj += alpha[i] * (G[i] + b[i]);

si.obj *= 0.5;

si.upper_bound_p = C[1];

si.upper_bound_n = C[0];

return true;

}

// return 1 if already optimal, return 0 otherwise

bool select_working_set( int& out_i, int& out_j )

{

// return i,j which maximize -grad(f)^T d , under constraint

// if alpha_i == C, d != +1

// if alpha_i == 0, d != -1

double Gmax1 = -DBL_MAX; // max { -grad(f)_i * d | y_i*d = +1 }

int Gmax1_idx = -1;

double Gmax2 = -DBL_MAX; // max { -grad(f)_i * d | y_i*d = -1 }

int Gmax2_idx = -1;

const schar* y = &y_vec[0];

const schar* alpha_status = &alpha_status_vec[0];

const double* G = &G_vec[0];

for( int i = 0; i < alpha_count; i++ )

{

double t;

if( y[i] > 0 ) // y = +1

{

if( !is_upper_bound(i) && (t = -G[i]) > Gmax1 ) // d = +1

{

Gmax1 = t;

Gmax1_idx = i;

}

if( !is_lower_bound(i) && (t = G[i]) > Gmax2 ) // d = -1

{

Gmax2 = t;

Gmax2_idx = i;

}

}

else // y = -1

{

if( !is_upper_bound(i) && (t = -G[i]) > Gmax2 ) // d = +1

{

Gmax2 = t;

Gmax2_idx = i;

}

if( !is_lower_bound(i) && (t = G[i]) > Gmax1 ) // d = -1

{

Gmax1 = t;

Gmax1_idx = i;

}

}

}

out_i = Gmax1_idx;

out_j = Gmax2_idx;

return Gmax1 + Gmax2 < eps;

}

void calc_rho( double& rho, double& r )

{

int nr_free = 0;

double ub = DBL_MAX, lb = -DBL_MAX, sum_free = 0;

const schar* y = &y_vec[0];

const schar* alpha_status = &alpha_status_vec[0];

const double* G = &G_vec[0];

for( int i = 0; i < alpha_count; i++ )

{

double yG = y[i]*G[i];

if( is_lower_bound(i) )

{

if( y[i] > 0 )

ub = MIN(ub,yG);

else

lb = MAX(lb,yG);

}

else if( is_upper_bound(i) )

{

if( y[i] < 0)

ub = MIN(ub,yG);

else

lb = MAX(lb,yG);

}

else

{

++nr_free;

sum_free += yG;

}

}

rho = nr_free > 0 ? sum_free/nr_free : (ub + lb)*0.5;

r = 0;

}

bool select_working_set_nu_svm( int& out_i, int& out_j )

{

// return i,j which maximize -grad(f)^T d , under constraint

// if alpha_i == C, d != +1

// if alpha_i == 0, d != -1

double Gmax1 = -DBL_MAX; // max { -grad(f)_i * d | y_i = +1, d = +1 }

int Gmax1_idx = -1;

double Gmax2 = -DBL_MAX; // max { -grad(f)_i * d | y_i = +1, d = -1 }

int Gmax2_idx = -1;

double Gmax3 = -DBL_MAX; // max { -grad(f)_i * d | y_i = -1, d = +1 }

int Gmax3_idx = -1;

double Gmax4 = -DBL_MAX; // max { -grad(f)_i * d | y_i = -1, d = -1 }

int Gmax4_idx = -1;

const schar* y = &y_vec[0];

const schar* alpha_status = &alpha_status_vec[0];

const double* G = &G_vec[0];

for( int i = 0; i < alpha_count; i++ )

{

double t;

if( y[i] > 0 ) // y == +1

{

if( !is_upper_bound(i) && (t = -G[i]) > Gmax1 ) // d = +1

{

Gmax1 = t;

Gmax1_idx = i;

}

if( !is_lower_bound(i) && (t = G[i]) > Gmax2 ) // d = -1

{

Gmax2 = t;

Gmax2_idx = i;

}

}

else // y == -1

{

if( !is_upper_bound(i) && (t = -G[i]) > Gmax3 ) // d = +1

{

Gmax3 = t;

Gmax3_idx = i;

}

if( !is_lower_bound(i) && (t = G[i]) > Gmax4 ) // d = -1

{

Gmax4 = t;

Gmax4_idx = i;

}

}

}

if( MAX(Gmax1 + Gmax2, Gmax3 + Gmax4) < eps )

return 1;

if( Gmax1 + Gmax2 > Gmax3 + Gmax4 )

{

out_i = Gmax1_idx;

out_j = Gmax2_idx;

}

else

{

out_i = Gmax3_idx;

out_j = Gmax4_idx;

}

return 0;

}

void calc_rho_nu_svm( double& rho, double& r )

{

int nr_free1 = 0, nr_free2 = 0;

double ub1 = DBL_MAX, ub2 = DBL_MAX;

double lb1 = -DBL_MAX, lb2 = -DBL_MAX;

double sum_free1 = 0, sum_free2 = 0;

const schar* y = &y_vec[0];

const schar* alpha_status = &alpha_status_vec[0];

const double* G = &G_vec[0];

for( int i = 0; i < alpha_count; i++ )

{

double G_i = G[i];

if( y[i] > 0 )

{

if( is_lower_bound(i) )

ub1 = MIN( ub1, G_i );

else if( is_upper_bound(i) )

lb1 = MAX( lb1, G_i );

else

{

++nr_free1;

sum_free1 += G_i;

}

}

else

{

if( is_lower_bound(i) )

ub2 = MIN( ub2, G_i );

else if( is_upper_bound(i) )

lb2 = MAX( lb2, G_i );

else

{

++nr_free2;

sum_free2 += G_i;