/*M///////////////////////////////////////////////////////////////////////////////////////

//

//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.

//

//  By downloading, copying, installing or using the software you agree to this license.

//  If you do not agree to this license, do not download, install,

//  copy or use the software.

//

//

//                        Intel License Agreement

//

// Copyright (C) 2000, Intel Corporation, all rights reserved.

// Third party copyrights are property of their respective owners.

//

// Redistribution and use in source and binary forms, with or without modification,

// are permitted provided that the following conditions are met:

//

//   * Redistribution's of source code must retain the above copyright notice,

//     this list of conditions and the following disclaimer.

//

//   * Redistribution's in binary form must reproduce the above copyright notice,

//     this list of conditions and the following disclaimer in the documentation

//     and/or other materials provided with the distribution.

//

//   * The name of Intel Corporation may not be used to endorse or promote products

//     derived from this software without specific prior written permission.

//

// This software is provided by the copyright holders and contributors "as is" and

// any express or implied warranties, including, but not limited to, the implied

// warranties of merchantability and fitness for a particular purpose are disclaimed.

// In no event shall the Intel Corporation or contributors be liable for any direct,

// indirect, incidental, special, exemplary, or consequential damages

// (including, but not limited to, procurement of substitute goods or services;

// loss of use, data, or profits; or business interruption) however caused

// and on any theory of liability, whether in contract, strict liability,

// or tort (including negligence or otherwise) arising in any way out of

// the use of this software, even if advised of the possibility of such damage.

//

//M*/

//Reformed function CvSVM::train_auto by Baruffaldi Juan Manuel -> [email protected]

#include "precomp.hpp"

/****************************************************************************************\

                                COPYRIGHT NOTICE

                                ----------------

  The code has been derived from libsvm library (version 2.6)

  (http://www.csie.ntu.edu.tw/~cjlin/libsvm).

  Here is the orignal copyright:

------------------------------------------------------------------------------------------

    Copyright (c) 2000-2003 Chih-Chung Chang and Chih-Jen Lin

    All rights reserved.

    Redistribution and use in source and binary forms, with or without

    modification, are permitted provided that the following conditions

    are met:

    1. Redistributions of source code must retain the above copyright

    notice, this list of conditions and the following disclaimer.

    2. Redistributions in binary form must reproduce the above copyright

    notice, this list of conditions and the following disclaimer in the

    documentation and/or other materials provided with the distribution.

    3. Neither name of copyright holders nor the names of its contributors

    may be used to endorse or promote products derived from this software

    without specific prior written permission.

    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS

    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT

    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR

    A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR

    CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,

    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,

    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR

    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF

    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING

    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS

    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

\****************************************************************************************/

using namespace cv;

#define CV_SVM_MIN_CACHE_SIZE  (40 << 20)  /* 40Mb */

#include <stdarg.h>

#include <ctype.h>

#if 1

typedef float Qfloat;

#define QFLOAT_TYPE CV_32F

#else

typedef double Qfloat;

#define QFLOAT_TYPE CV_64F

#endif

// Param Grid

bool CvParamGrid::check() const

{

    bool ok = false;

    CV_FUNCNAME( "CvParamGrid::check" );

    __BEGIN__;

    if( min_val > max_val )

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

    if( min_val < DBL_EPSILON )

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

    if( step < 1. + FLT_EPSILON )

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

    ok = true;

    __END__;

    return ok;

}

CvParamGrid CvSVM::get_default_grid( int param_id )

{

    CvParamGrid grid;

    if( param_id == CvSVM::C )

    {

        grid.min_val = 0.1;

        grid.max_val = 500;

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

    }

    else if( param_id == CvSVM::GAMMA )

    {

        grid.min_val = 1e-5;

        grid.max_val = 0.6;

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

    }

    else if( param_id == CvSVM::P )

    {

        grid.min_val = 0.01;

        grid.max_val = 100;

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

    }

    else if( param_id == CvSVM::NU )

    {

        grid.min_val = 0.01;

        grid.max_val = 0.2;

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

    }

    else if( param_id == CvSVM::COEF )

    {

        grid.min_val = 0.1;

        grid.max_val = 300;

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

    }

    else if( param_id == CvSVM::DEGREE )

    {

        grid.min_val = 0.01;

        grid.max_val = 4;

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

    }

    else

        cvError( CV_StsBadArg, "CvSVM::get_default_grid", "Invalid type of parameter "

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

    return grid;

}

// SVM training parameters

CvSVMParams::CvSVMParams() :

    svm_type(CvSVM::C_SVC), kernel_type(CvSVM::RBF), degree(0),

    gamma(1), coef0(0), C(1), nu(0), p(0), class_weights(0)

{

    term_crit = cvTermCriteria( CV_TERMCRIT_ITER+CV_TERMCRIT_EPS, 1000, FLT_EPSILON );

}

CvSVMParams::CvSVMParams( int _svm_type, int _kernel_type,

    double _degree, double _gamma, double _coef0,

    double _Con, double _nu, double _p,

    CvMat* _class_weights, CvTermCriteria _term_crit ) :

    svm_type(_svm_type), kernel_type(_kernel_type),

    degree(_degree), gamma(_gamma), coef0(_coef0),

    C(_Con), nu(_nu), p(_p), class_weights(_class_weights), term_crit(_term_crit)

{

}

/////////////////////////////////////// SVM kernel ///////////////////////////////////////

CvSVMKernel::CvSVMKernel()

{

    clear();

}

void CvSVMKernel::clear()

{

    params = 0;

    calc_func = 0;

}

CvSVMKernel::~CvSVMKernel()

{

}

CvSVMKernel::CvSVMKernel( const CvSVMParams* _params, Calc _calc_func )

{

    clear();

    create( _params, _calc_func );

}

bool CvSVMKernel::create( const CvSVMParams* _params, Calc _calc_func )

{

    clear();

    params = _params;

    calc_func = _calc_func;

    if( !calc_func )

        calc_func = params->kernel_type == CvSVM::RBF ? &CvSVMKernel::calc_rbf :

                    params->kernel_type == CvSVM::POLY ? &CvSVMKernel::calc_poly :

                    params->kernel_type == CvSVM::SIGMOID ? &CvSVMKernel::calc_sigmoid :

                    &CvSVMKernel::calc_linear;

    return true;

}

void CvSVMKernel::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];

        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 CvSVMKernel::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 CvSVMKernel::calc_poly( int vcount, int var_count, const float** vecs,

                             const float* another, Qfloat* results )

{

    CvMat R = cvMat( 1, vcount, QFLOAT_TYPE, results );

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

    if( vcount > 0 )

        cvPow( &R, &R, params->degree );

}

void CvSVMKernel::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];

        double e = exp(-fabs(t));

        if( t > 0 )

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

        else

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

    }

}

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

                            const float* another, Qfloat* results )

{

    CvMat R = cvMat( 1, vcount, QFLOAT_TYPE, results );

    double gamma = -params->gamma;

    int j, k;

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

    {

        const float* sample = vecs[j];

        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 )

        cvExp( &R, &R );

}

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

                        const float* another, Qfloat* results )

{

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

    int j;

    (this->*calc_func)( vcount, var_count, vecs, another, results );

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

    {

        if( results[j] > max_val )

            results[j] = max_val;

    }

}

// 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

//

void CvSVMSolver::clear()

{

    G = 0;

    alpha = 0;

    y = 0;

    b = 0;

    buf[0] = buf[1] = 0;

    cvReleaseMemStorage( &storage );

    kernel = 0;

    select_working_set_func = 0;

    calc_rho_func = 0;

    rows = 0;

    samples = 0;

    get_row_func = 0;

}

CvSVMSolver::CvSVMSolver()

{

    storage = 0;

    clear();

}

CvSVMSolver::~CvSVMSolver()

{

    clear();

}

CvSVMSolver::CvSVMSolver( int _sample_count, int _var_count, const float** _samples, schar* _y,

                int _alpha_count, double* _alpha, double _Cp, double _Cn,

                CvMemStorage* _storage, CvSVMKernel* _kernel, GetRow _get_row,

                SelectWorkingSet _select_working_set, CalcRho _calc_rho )

{

    storage = 0;

    create( _sample_count, _var_count, _samples, _y, _alpha_count, _alpha, _Cp, _Cn,

            _storage, _kernel, _get_row, _select_working_set, _calc_rho );

}

bool CvSVMSolver::create( int _sample_count, int _var_count, const float** _samples, schar* _y,

                int _alpha_count, double* _alpha, double _Cp, double _Cn,

                CvMemStorage* _storage, CvSVMKernel* _kernel, GetRow _get_row,

                SelectWorkingSet _select_working_set, CalcRho _calc_rho )

{

    bool ok = false;

    int i, svm_type;

    CV_FUNCNAME( "CvSVMSolver::create" );

    __BEGIN__;

    int rows_hdr_size;

    clear();

    sample_count = _sample_count;

    var_count = _var_count;

    samples = _samples;

    y = _y;

    alpha_count = _alpha_count;

    alpha = _alpha;

    kernel = _kernel;

    C[0] = _Cn;

    C[1] = _Cp;

    eps = kernel->params->term_crit.epsilon;

    max_iter = kernel->params->term_crit.max_iter;

    storage = cvCreateChildMemStorage( _storage );

    b = (double*)cvMemStorageAlloc( storage, alpha_count*sizeof(b[0]));

    alpha_status = (schar*)cvMemStorageAlloc( storage, alpha_count*sizeof(alpha_status[0]));

    G = (double*)cvMemStorageAlloc( storage, alpha_count*sizeof(G[0]));

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

        buf[i] = (Qfloat*)cvMemStorageAlloc( storage, sample_count*2*sizeof(buf[i][0]) );

    svm_type = kernel->params->svm_type;

    select_working_set_func = _select_working_set;

    if( !select_working_set_func )

        select_working_set_func = svm_type == CvSVM::NU_SVC || svm_type == CvSVM::NU_SVR ?

        &CvSVMSolver::select_working_set_nu_svm : &CvSVMSolver::select_working_set;

    calc_rho_func = _calc_rho;

    if( !calc_rho_func )

        calc_rho_func = svm_type == CvSVM::NU_SVC || svm_type == CvSVM::NU_SVR ?

            &CvSVMSolver::calc_rho_nu_svm : &CvSVMSolver::calc_rho;

    get_row_func = _get_row;

    if( !get_row_func )

        get_row_func = params->svm_type == CvSVM::EPS_SVR ||

                       params->svm_type == CvSVM::NU_SVR ? &CvSVMSolver::get_row_svr :

                       params->svm_type == CvSVM::C_SVC ||

                       params->svm_type == CvSVM::NU_SVC ? &CvSVMSolver::get_row_svc :

                       &CvSVMSolver::get_row_one_class;

    cache_line_size = sample_count*sizeof(Qfloat);

    // cache size = max(num_of_samples^2*sizeof(Qfloat)*0.25, 64Kb)

    // (assuming that for large training sets ~25% of Q matrix is used)

    cache_size = MAX( cache_line_size*sample_count/4, CV_SVM_MIN_CACHE_SIZE );

    // the size of Q matrix row headers

    rows_hdr_size = sample_count*sizeof(rows[0]);

    if( rows_hdr_size > storage->block_size )

        CV_ERROR( CV_StsOutOfRange, "Too small storage block size" );

    lru_list.prev = lru_list.next = &lru_list;

    rows = (CvSVMKernelRow*)cvMemStorageAlloc( storage, rows_hdr_size );

    memset( rows, 0, rows_hdr_size );

    ok = true;

    __END__;

    return ok;

}

float* CvSVMSolver::get_row_base( int i, bool* _existed )

{

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

    CvSVMKernelRow* row = rows + i1;

    bool existed = row->data != 0;

    Qfloat* data;

    if( existed || cache_size <= 0 )

    {

        CvSVMKernelRow* del_row = existed ? row : lru_list.prev;

        data = del_row->data;

        assert( data != 0 );

        // delete row from the LRU list

        del_row->data = 0;

        del_row->prev->next = del_row->next;

        del_row->next->prev = del_row->prev;

    }

    else

    {

        data = (Qfloat*)cvMemStorageAlloc( storage, cache_line_size );

        cache_size -= cache_line_size;

    }

    // insert row into the LRU list

    row->data = data;

    row->prev = &lru_list;

    row->next = lru_list.next;

    row->prev->next = row->next->prev = row;

    if( !existed )

    {

        kernel->calc( sample_count, var_count, samples, samples[i1], row->data );

    }

    if( _existed )

        *_existed = existed;

    return row->data;

}

float* CvSVMSolver::get_row_svc( int i, float* row, float*, bool existed )

{

    if( !existed )

    {

        const schar* _y = y;

        int j, len = sample_count;

        assert( _y && i < 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;

}

float* CvSVMSolver::get_row_one_class( int, float* row, float*, bool )

{

    return row;

}

float* CvSVMSolver::get_row_svr( int i, float* row, float* dst, bool )

{

    int j, len = sample_count;

    Qfloat* dst_pos = dst;

    Qfloat* dst_neg = dst + len;

    if( i >= len )

    {

        Qfloat* temp;

        CV_SWAP( dst_pos, dst_neg, temp );

    }

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

    {

        Qfloat t = row[j];

        dst_pos[j] = t;

        dst_neg[j] = -t;

    }

    return dst;

}

float* CvSVMSolver::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 CvSVMSolver::solve_generic( CvSVMSolutionInfo& si )

{

    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] );

            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] );

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

        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

CvSVMSolver::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;

    int i;

    for( 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

CvSVMSolver::calc_rho( double& rho, double& r )

{

    int i, nr_free = 0;

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

    for( 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

CvSVMSolver::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;

    int i;

    for( 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

CvSVMSolver::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;

    double r1, r2;

    int i;

    for( 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;

            }

        }

    }

    r1 = nr_free1 > 0 ? sum_free1/nr_free1 : (ub1 + lb1)*0.5;

    r2 = nr_free2 > 0 ? sum_free2/nr_free2 : (ub2 + lb2)*0.5;

    rho = (r1 - r2)*0.5;

    r = (r1 + r2)*0.5;

}

/*

///////////////////////// construct and solve various formulations ///////////////////////

*/

bool CvSVMSolver::solve_c_svc( int _sample_count, int _var_count, const float** _samples, schar* _y,

                               double _Cp, double _Cn, CvMemStorage* _storage,

                               CvSVMKernel* _kernel, double* _alpha, CvSVMSolutionInfo& _si )

{

    int i;

    if( !create( _sample_count, _var_count, _samples, _y, _sample_count,

                 _alpha, _Cp, _Cn, _storage, _kernel, &CvSVMSolver::get_row_svc,

                 &CvSVMSolver::select_working_set, &CvSVMSolver::calc_rho ))

        return false;

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

    {

        alpha[i] = 0;

        b[i] = -1;

    }

    if( !solve_generic( _si ))

        return false;

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

        alpha[i] *= y[i];

    return true;

}

bool CvSVMSolver::solve_nu_svc( int _sample_count, int _var_count, const float** _samples, schar* _y,

                                CvMemStorage* _storage, CvSVMKernel* _kernel,

                                double* _alpha, CvSVMSolutionInfo& _si )

{

    int i;

    double sum_pos, sum_neg, inv_r;

    if( !create( _sample_count, _var_count, _samples, _y, _sample_count,

                 _alpha, 1., 1., _storage, _kernel, &CvSVMSolver::get_row_svc,

                 &CvSVMSolver::select_working_set_nu_svm, &CvSVMSolver::calc_rho_nu_svm ))

        return false;

    sum_pos = kernel->params->nu * sample_count * 0.5;

    sum_neg = kernel->params->nu * sample_count * 0.5;

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

    {

        if( y[i] > 0 )

        {

            alpha[i] = MIN(1.0, sum_pos);

            sum_pos -= alpha[i];

        }

        else

        {

            alpha[i] = MIN(1.0, sum_neg);

            sum_neg -= alpha[i];

        }

        b[i] = 0;

    }

    if( !solve_generic( _si ))

        return false;

    inv_r = 1./_si.r;

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

        alpha[i] *= y[i]*inv_r;

    _si.rho *= inv_r;

    _si.obj *= (inv_r*inv_r);

    _si.upper_bound_p = inv_r;

    _si.upper_bound_n = inv_r;

    return true;

}

bool CvSVMSolver::solve_one_class( int _sample_count, int _var_count, const float** _samples,

                                   CvMemStorage* _storage, CvSVMKernel* _kernel,

                                   double* _alpha, CvSVMSolutionInfo& _si )

{

    int i, n;

    double nu = _kernel->params->nu;

    if( !create( _sample_count, _var_count, _samples, 0, _sample_count,

                 _alpha, 1., 1., _storage, _kernel, &CvSVMSolver::get_row_one_class,

                 &CvSVMSolver::select_working_set, &CvSVMSolver::calc_rho ))

        return false;

    y = (schar*)cvMemStorageAlloc( storage, sample_count*sizeof(y[0]) );

    n = cvRound( nu*sample_count );

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

    {

        y[i] = 1;

        b[i] = 0;

        alpha[i] = i < n ? 1 : 0;

    }

    if( n < sample_count )

        alpha[n] = nu * sample_count - n;

    else

        alpha[n-1] = nu * sample_count - (n-1);

    return solve_generic(_si);

}

bool CvSVMSolver::solve_eps_svr( int _sample_count, int _var_count, const float** _samples,

                                 const float* _y, CvMemStorage* _storage,

                                 CvSVMKernel* _kernel, double* _alpha, CvSVMSolutionInfo& _si )

{

    int i;

    double p = _kernel->params->p, kernel_param_c = _kernel->params->C;

    if( !create( _sample_count, _var_count, _samples, 0,

                 _sample_count*2, 0, kernel_param_c, kernel_param_c, _storage, _kernel, &CvSVMSolver::get_row_svr,

                 &CvSVMSolver::select_working_set, &CvSVMSolver::calc_rho ))

        return false;

    y = (schar*)cvMemStorageAlloc( storage, sample_count*2*sizeof(y[0]) );

    alpha = (double*)cvMemStorageAlloc( storage, alpha_count*sizeof(alpha[0]) );

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

    {

        alpha[i] = 0;

        b[i] = p - _y[i];

        y[i] = 1;

        alpha[i+sample_count] = 0;

        b[i+sample_count] = p + _y[i];

        y[i+sample_count] = -1;

    }

    if( !solve_generic( _si ))

        return false;

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

        _alpha[i] = alpha[i] - alpha[i+sample_count];

    return true;

}

bool CvSVMSolver::solve_nu_svr( int _sample_count, int _var_count, const float** _samples,

                                const float* _y, CvMemStorage* _storage,

                                CvSVMKernel* _kernel, double* _alpha, CvSVMSolutionInfo& _si )

{

    int i;

    double kernel_param_c = _kernel->params->C, sum;

    if( !create( _sample_count, _var_count, _samples, 0,

                 _sample_count*2, 0, 1., 1., _storage, _kernel, &CvSVMSolver::get_row_svr,

                 &CvSVMSolver::select_working_set_nu_svm, &CvSVMSolver::calc_rho_nu_svm ))

        return false;

    y = (schar*)cvMemStorageAlloc( storage, sample_count*2*sizeof(y[0]) );

    alpha = (double*)cvMemStorageAlloc( storage, alpha_count*sizeof(alpha[0]) );

    sum = kernel_param_c * _kernel->params->nu * sample_count * 0.5;

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

    {

        alpha[i] = alpha[i + sample_count] = MIN(sum, kernel_param_c);

        sum -= alpha[i];

        b[i] = -_y[i];

        y[i] = 1;

        b[i + sample_count] = _y[i];

        y[i + sample_count] = -1;

    }

    if( !solve_generic( _si ))

        return false;

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

        _alpha[i] = alpha[i] - alpha[i+sample_count];

    return true;

}

//////////////////////////////////////////////////////////////////////////////////////////

CvSVM::CvSVM()

{

    decision_func = 0;

    class_labels = 0;

    class_weights = 0;

    storage = 0;

    var_idx = 0;

    kernel = 0;

    solver = 0;

    default_model_name = "my_svm";

    clear();

}

CvSVM::~CvSVM()

{

    clear();

}

void CvSVM::clear()

{

    cvFree( &decision_func );

    cvReleaseMat( &class_labels );

    cvReleaseMat( &class_weights );

    cvReleaseMemStorage( &storage );

    cvReleaseMat( &var_idx );

    delete kernel;

    delete solver;

    kernel = 0;

    solver = 0;

    var_all = 0;

    sv = 0;

    sv_total = 0;

}

CvSVM::CvSVM( const CvMat* _train_data, const CvMat* _responses,

    const CvMat* _var_idx, const CvMat* _sample_idx, CvSVMParams _params )

{

    decision_func = 0;

    class_labels = 0;

    class_weights = 0;

    storage = 0;

    var_idx = 0;

    kernel = 0;

    solver = 0;

    default_model_name = "my_svm";

    train( _train_data, _responses, _var_idx, _sample_idx, _params );

}

int CvSVM::get_support_vector_count() const

{

    return sv_total;

}

const float* CvSVM::get_support_vector(int i) const

{

    return sv && (unsigned)i < (unsigned)sv_total ? sv[i] : 0;

}

bool CvSVM::set_params( const CvSVMParams& _params )

{

    bool ok = false;

    CV_FUNCNAME( "CvSVM::set_params" );

    __BEGIN__;

    int kernel_type, svm_type;

    params = _params;

    kernel_type = params.kernel_type;

    svm_type = params.svm_type;

    if( kernel_type != LINEAR && kernel_type != POLY &&

        kernel_type != SIGMOID && kernel_type != RBF )

        CV_ERROR( CV_StsBadArg, "Unknown/unsupported kernel type" );

    if( kernel_type == LINEAR )

        params.gamma = 1;

    else if( params.gamma <= 0 )

        CV_ERROR( CV_StsOutOfRange, "gamma parameter of the kernel must be positive" );

    if( kernel_type != SIGMOID && kernel_type != POLY )

        params.coef0 = 0;

    else if( params.coef0 < 0 )

        CV_ERROR( CV_StsOutOfRange, "The kernel parameter <coef0> must be positive or zero" );

    if( kernel_type != POLY )

        params.degree = 0;

    else if( params.degree <= 0 )

        CV_ERROR( CV_StsOutOfRange, "The kernel parameter <degree> must be positive" );

    if( svm_type != C_SVC && svm_type != NU_SVC &&

        svm_type != ONE_CLASS && svm_type != EPS_SVR &&

        svm_type != NU_SVR )

        CV_ERROR( CV_StsBadArg, "Unknown/unsupported SVM type" );

    if( svm_type == ONE_CLASS || svm_type == NU_SVC )

        params.C = 0;

    else if( params.C <= 0 )

        CV_ERROR( CV_StsOutOfRange, "The parameter C must be positive" );

    if( svm_type == C_SVC || svm_type == EPS_SVR )

        params.nu = 0;

    else if( params.nu <= 0 || params.nu >= 1 )

        CV_ERROR( CV_StsOutOfRange, "The parameter nu must be between 0 and 1" );

    if( svm_type != EPS_SVR )

        params.p = 0;

    else if( params.p <= 0 )

        CV_ERROR( CV_StsOutOfRange, "The parameter p must be positive" );

    if( svm_type != C_SVC )

        params.class_weights = 0;

    params.term_crit = cvCheckTermCriteria( params.term_crit, DBL_EPSILON, INT_MAX );

    params.term_crit.epsilon = MAX( params.term_crit.epsilon, DBL_EPSILON );

    ok = true;

    __END__;

    return ok;

}

void CvSVM::create_kernel()

{

    kernel = new CvSVMKernel(&params,0);

}

void CvSVM::create_solver( )

{

    solver = new CvSVMSolver;

}

// switching function

bool CvSVM::train1( int sample_count, int var_count, const float** samples,

                    const void* _responses, double Cp, double Cn,

                    CvMemStorage* _storage, double* alpha, double& rho )

{

    bool ok = false;

    //CV_FUNCNAME( "CvSVM::train1" );

    __BEGIN__;

    CvSVMSolutionInfo si;

    int svm_type = params.svm_type;

    si.rho = 0;

    ok = svm_type == C_SVC ? solver->solve_c_svc( sample_count, var_count, samples, (schar*)_responses,

                                                  Cp, Cn, _storage, kernel, alpha, si ) :

         svm_type == NU_SVC ? solver->solve_nu_svc( sample_count, var_count, samples, (schar*)_responses,

                                                    _storage, kernel, alpha, si ) :

         svm_type == ONE_CLASS ? solver->solve_one_class( sample_count, var_count, samples,

                                                          _storage, kernel, alpha, si ) :

         svm_type == EPS_SVR ? solver->solve_eps_svr( sample_count, var_count, samples, (float*)_responses,

                                                      _storage, kernel, alpha, si ) :

         svm_type == NU_SVR ? solver->solve_nu_svr( sample_count, var_count, samples, (float*)_responses,

                                                    _storage, kernel, alpha, si ) : false;

    rho = si.rho;

    __END__;

    return ok;

}

bool CvSVM::do_train( int svm_type, int sample_count, int var_count, const float** samples,

                    const CvMat* responses, CvMemStorage* temp_storage, double* alpha )

{

    bool ok = false;

    CV_FUNCNAME( "CvSVM::do_train" );

    __BEGIN__;

    CvSVMDecisionFunc* df = 0;

    const int sample_size = var_count*sizeof(samples[0][0]);

    int i, j, k;

    cvClearMemStorage( storage );

    if( svm_type == ONE_CLASS || svm_type == EPS_SVR || svm_type == NU_SVR )

    {

        int sv_count = 0;

        CV_CALL( decision_func = df =

            (CvSVMDecisionFunc*)cvAlloc( sizeof(df[0]) ));

        df->rho = 0;

        if( !train1( sample_count, var_count, samples, svm_type == ONE_CLASS ? 0 :

            responses->data.i, 0, 0, temp_storage, alpha, df->rho ))

            EXIT;

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

            sv_count += fabs(alpha[i]) > 0;

        CV_Assert(sv_count != 0);

        sv_total = df->sv_count = sv_count;

        CV_CALL( df->alpha = (double*)cvMemStorageAlloc( storage, sv_count*sizeof(df->alpha[0])) );

        CV_CALL( sv = (float**)cvMemStorageAlloc( storage, sv_count*sizeof(sv[0])));

        for( i = k = 0; i < sample_count; i++ )

        {

            if( fabs(alpha[i]) > 0 )

            {

                CV_CALL( sv[k] = (float*)cvMemStorageAlloc( storage, sample_size ));

                memcpy( sv[k], samples[i], sample_size );

                df->alpha[k++] = alpha[i];

            }

        }

    }

    else

    {

        int class_count = class_labels->cols;

        int* sv_tab = 0;

        const float** temp_samples = 0;

        int* class_ranges = 0;