各位大大好,小弟從web抓到一個C++的SVM分類程式，但卻無法進行combination成EXE檔，程式碼如下，請前輩們參照，謝謝您們  

 
#include 
#include 
#include 
#include 
#include "svm.h"
#define Malloc(type,n) (type *)malloc((n)*sizeof(type))    void exit_with_help()
{
        printf(
        "Usage: svm-train [options] training_set_file [model_file]\n"
        "options:\n"
        "-s svm_type : set type of SVM (default 0)\n"
        "        0 -- C-SVC\n"
        "        1 -- nu-SVC\n"
        "        2 -- one-class SVM\n"
        "        3 -- epsilon-SVR\n"
        "        4 -- nu-SVR\n"
        "-t kernel_type : set type of kernel function (default 2)\n"
        "        0 -- linear: u'*v\n"
        "        1 -- polynomial: (gamma*u'*v   coef0)^degree\n"
        "        2 -- radial basis function: exp(-gamma*|u-v|^2)\n"
        "        3 -- sigmoid: tanh(gamma*u'*v   coef0)\n"
        "-d degree : set degree in kernel function (default 3)\n"
        "-g gamma : set gamma in kernel function (default 1/k)\n"
        "-r coef0 : set coef0 in kernel function (default 0)\n"
        "-c cost : set the parameter C of C-SVC, epsilon-SVR, and nu-SVR (default 1)\n"
        "-n nu : set the parameter nu of nu-SVC, one-class SVM, and nu-SVR (default 0.5)\n"
        "-p epsilon : set the epsilon in loss function of epsilon-SVR (default 0.1)\n"
        "-m cachesize : set cache memory size in MB (default 40)\n"
        "-e epsilon : set tolerance of termination criterion (default 0.001)\n"
        "-h shrinking: whether to use the shrinking heuristics, 0 or 1 (default 1)\n"
        "-b probability_estimates: whether to train a SVC or SVR model for probability estimates, 0 or 1 (default 0)\n"
        "-wi weight: set the parameter C of class i to weight*C, for C-SVC (default 1)\n"
        "-v n: n-fold cross validation mode\n"
        );
        exit(1);
}    void parse_command_line(int argc, char **argv, char *input_file_name, char *model_file_name);
void read_problem(const char *filename);
void do_cross_validation();    struct svm_parameter param;                // set by parse_command_line
struct svm_problem prob;                // set by read_problem
struct svm_model *model;
struct svm_node *x_space;
int cross_validation;
int nr_fold;    int main(int argc, char **argv)
{
        char input_file_name[1024];
        char model_file_name[1024];
        const char *error_msg;            parse_command_line(argc, argv, input_file_name, model_file_name);
        read_problem(input_file_name);
        error_msg = svm_check_parameter(&prob,¶m);            if(error_msg)
        {
                fprintf(stderr,"Error: %s\n",error_msg);
                exit(1);
        }            if(cross_validation)
        {
                do_cross_validation();
        }
        else
        {
                model = svm_train(&prob,¶m);
                svm_save_model(model_file_name,model);
                svm_destroy_model(model);
        }
        svm_destroy_param(¶m);
        free(prob.y);
        free(prob.x);
        free(x_space);            return 0;
}    void do_cross_validation()
{
        int i;
        int total_correct = 0;
        double total_error = 0;
        double sumv = 0, sumy = 0, sumvv = 0, sumyy = 0, sumvy = 0;
        double *target = Malloc(double,prob.l);            svm_cross_validation(&prob,¶m,nr_fold,target);
        if(param.svm_type == EPSILON_SVR ||
           param.svm_type == NU_SVR)
        {
                for(i=0;i=argc)
                        exit_with_help();
                switch(argv[i-1][1])
                {
                        case 's':
                                param.svm_type = atoi(argv[i]);
                                break;
                        case 't':
                                param.kernel_type = atoi(argv[i]);
                                break;
                        case 'd':
                                param.degree = atof(argv[i]);
                                break;
                        case 'g':
                                param.gamma = atof(argv[i]);
                                break;
                        case 'r':
                                param.coef0 = atof(argv[i]);
                                break;
                        case 'n':
                                param.nu = atof(argv[i]);
                                break;
                        case 'm':
                                param.cache_size = atof(argv[i]);
                                break;
                        case 'c':
                                param.C = atof(argv[i]);
                                break;
                        case 'e':
                                param.eps = atof(argv[i]);
                                break;
                        case 'p':
                                param.p = atof(argv[i]);
                                break;
                        case 'h':
                                param.shrinking = atoi(argv[i]);
                                break;
                        case 'b':
                                param.probability = atoi(argv[i]);
                                break;
                        case 'v':
                                cross_validation = 1;
                                nr_fold = atoi(argv[i]);
                                if(nr_fold < 2)
                                {
                                        fprintf(stderr,"n-fold cross validation: n must >= 2\n");
                                        exit_with_help();
                                }
                                break;
                        case 'w':
                                  param.nr_weight;
                                param.weight_label = (int *)realloc(param.weight_label,sizeof(int)*param.nr_weight);
                                param.weight = (double *)realloc(param.weight,sizeof(double)*param.nr_weight);
                                param.weight_label[param.nr_weight-1] = atoi(&argv[i-1][2]);
                                param.weight[param.nr_weight-1] = atof(argv[i]);
                                break;
                        default:
                                fprintf(stderr,"unknown option\n");
                                exit_with_help();
                }
        }            // determine filenames            if(i>=argc)
                exit_with_help();            strcpy(input_file_name, argv[i]);            if(i max_index)
                        max_index = x_space[j-1].index;
                x_space[j  ].index = -1;
        }            if(param.gamma == 0)
                param.gamma = 1.0/max_index;            fclose(fp);
}

發表人 - oooopppp1 於 2005/06/20 15:30:05

您好:        PO程式碼的方式請參考版規說明,煩請修改謝謝您的配合 >

where is your header file "svm.h", and how could we debug it ? -----------------------------------------------

Compile你提供的程式碼大致有兩個問題！ 1.你的程式碼不完整，例如在svm.h裏宣告的很多函數(非內建函數)，在主程式裡並找不到它的定義，因此呼叫時產生錯誤，例如svm_...等開頭的函數通通都沒有！ 2.程式碼裡有不可見字元，例如   error_msg = svm_check_parameter(&prob,¶m); // compiler不知道如何解譯！ 不可見字元在很多函數裡都有喔！ 因此推斷你download到的只是一個不完整的範例。 ----------------------------------------------- Creation is the fundation of promotion. 發表人 - stallion 於 2005/06/21 20:49:46

首先謝謝各位大大的回應，一切錯誤來源來自弟沒把另一個code併入，結論的確是如大大所述並不完整，需要將另一段副程式add加入專案如此才算完整，弟在此將今天剛找到的另一程式碼附上，並將此問題結案，再次謝謝大大門的鼎力相助，謝謝您們^_^,the code as follow.. ---  

 
#include <math.h>
#include 
#include 
#include 
#include 
#include 
#include 
#include "svm.h"
typedef float Qfloat;
typedef signed char schar;
#ifndef min
template  inline T min(T x,T y) { return (x inline T max(T x,T y) { return (x>y)?x:y; }
#endif
template  inline void swap(T& x, T& y) { T t=x; x=y; y=t; }
template  inline void clone(T*& dst, S* src, int n)
{
        dst = new T[n];
        memcpy((void *)dst,(void *)src,sizeof(T)*n);
}
#define INF HUGE_VAL
#define TAU 1e-12
#define Malloc(type,n) (type *)malloc((n)*sizeof(type))
#if 1
void info(char *fmt,...)
{
        va_list ap;
        va_start(ap,fmt);
        vprintf(fmt,ap);
        va_end(ap);
}
void info_flush()
{
        fflush(stdout);
}
#else
void info(char *fmt,...) {}
void info_flush() {}
#endif    //
// Kernel Cache
//
// l is the number of total data items
// size is the cache size limit in bytes
//
class Cache
{
public:
        Cache(int l,int size);
        ~Cache();            // request data [0,len)
        // return some position p where [p,len) need to be filled
        // (p >= len if nothing needs to be filled)
        int get_data(const int index, Qfloat **data, int len);
        void swap_index(int i, int j);        // future_option
private:
        int l;
        int size;
        struct head_t
        {
                head_t *prev, *next;        // a cicular list
                Qfloat *data;
                int len;                // data[0,len) is cached in this entry
        };            head_t *head;
        head_t lru_head;
        void lru_delete(head_t *h);
        void lru_insert(head_t *h);
};    Cache::Cache(int l_,int size_):l(l_),size(size_)
{
        head = (head_t *)calloc(l,sizeof(head_t));        // initialized to 0
        size /= sizeof(Qfloat);
        size -= l * sizeof(head_t) / sizeof(Qfloat);
        size = max(size, 2*l);        // cache must be large enough for two columns
        lru_head.next = lru_head.prev = &lru_head;
}    Cache::~Cache()
{
        for(head_t *h = lru_head.next; h != &lru_head; h=h->next)
                free(h->data);
        free(head);
}    void Cache::lru_delete(head_t *h)
{
        // delete from current location
        h->prev->next = h->next;
        h->next->prev = h->prev;
}    void Cache::lru_insert(head_t *h)
{
        // insert to last position
        h->next = &lru_head;
        h->prev = lru_head.prev;
        h->prev->next = h;
        h->next->prev = h;
}    int Cache::get_data(const int index, Qfloat **data, int len)
{
        head_t *h = &head[index];
        if(h->len) lru_delete(h);
        int more = len - h->len;            if(more > 0)
        {
                // free old space
                while(size < more)
                {
                        head_t *old = lru_head.next;
                        lru_delete(old);
                        free(old->data);
                        size  = old->len;
                        old->data = 0;
                        old->len = 0;
                }                    // allocate new space
                h->data = (Qfloat *)realloc(h->data,sizeof(Qfloat)*len);
                size -= more;
                swap(h->len,len);
        }            lru_insert(h);
        *data = h->data;
        return len;
}    void Cache::swap_index(int i, int j)
{
        if(i==j) return;            if(head[i].len) lru_delete(&head[i]);
        if(head[j].len) lru_delete(&head[j]);
        swap(head[i].data,head[j].data);
        swap(head[i].len,head[j].len);
        if(head[i].len) lru_insert(&head[i]);
        if(head[j].len) lru_insert(&head[j]);            if(i>j) swap(i,j);
        for(head_t *h = lru_head.next; h!=&lru_head; h=h->next)
        {
                if(h->len > i)
                {
                        if(h->len > j)
                                swap(h->data[i],h->data[j]);
                        else
                        {
                                // give up
                                lru_delete(h);
                                free(h->data);
                                size  = h->len;
                                h->data = 0;
                                h->len = 0;
                        }
                }
        }
}    //
// Kernel evaluation
//
// the static method k_function is for doing single kernel evaluation
// the constructor of Kernel prepares to calculate the l*l kernel matrix
// the member function get_Q is for getting one column from the Q Matrix
//
class QMatrix {
public:
        virtual Qfloat *get_Q(int column, int len) const = 0;
        virtual Qfloat *get_QD() const = 0;
        virtual void swap_index(int i, int j) const = 0;
};    class Kernel: public QMatrix {
public:
        Kernel(int l, svm_node * const * x, const svm_parameter& param);
        virtual ~Kernel();            static double k_function(const svm_node *x, const svm_node *y,
                                 const svm_parameter& param);
        virtual Qfloat *get_Q(int column, int len) const = 0;
        virtual Qfloat *get_QD() const = 0;
        virtual void swap_index(int i, int j) const        // no so const...
        {
                swap(x[i],x[j]);
                if(x_square) swap(x_square[i],x_square[j]);
        }
protected:            double (Kernel::*kernel_function)(int i, int j) const;    private:
        const svm_node **x;
        double *x_square;            // svm_parameter
        const int kernel_type;
        const double degree;
        const double gamma;
        const double coef0;            static double dot(const svm_node *px, const svm_node *py);
        double kernel_linear(int i, int j) const
        {
                return dot(x[i],x[j]);
        }
        double kernel_poly(int i, int j) const
        {
                return pow(gamma*dot(x[i],x[j]) coef0,degree);
        }
        double kernel_rbf(int i, int j) const
        {
                return exp(-gamma*(x_square[i] x_square[j]-2*dot(x[i],x[j])));
        }
        double kernel_sigmoid(int i, int j) const
        {
                return tanh(gamma*dot(x[i],x[j]) coef0);
        }
};    Kernel::Kernel(int l, svm_node * const * x_, const svm_parameter& param)
:kernel_type(param.kernel_type), degree(param.degree),
 gamma(param.gamma), coef0(param.coef0)
{
        switch(kernel_type)
        {
                case LINEAR:
                        kernel_function = &Kernel::kernel_linear;
                        break;
                case POLY:
                        kernel_function = &Kernel::kernel_poly;
                        break;
                case RBF:
                        kernel_function = &Kernel::kernel_rbf;
                        break;
                case SIGMOID:
                        kernel_function = &Kernel::kernel_sigmoid;
                        break;
        }            clone(x,x_,l);            if(kernel_type == RBF)
        {
                x_square = new double[l];
                for(int i=0;iindex != -1 && py->index != -1)
        {
                if(px->index == py->index)
                {
                        sum  = px->value * py->value;
                          px;
                          py;
                }
                else
                {
                        if(px->index > py->index)
                                  py;
                        else
                                  px;
                }                        
        }
        return sum;
}    double Kernel::k_function(const svm_node *x, const svm_node *y,
                          const svm_parameter& param)
{
        switch(param.kernel_type)
        {
                case LINEAR:
                        return dot(x,y);
                case POLY:
                        return pow(param.gamma*dot(x,y) param.coef0,param.degree);
                case RBF:
                {
                        double sum = 0;
                        while(x->index != -1 && y->index !=-1)
                        {
                                if(x->index == y->index)
                                {
                                        double d = x->value - y->value;
                                        sum  = d*d;
                                          x;
                                          y;
                                }
                                else
                                {
                                        if(x->index > y->index)
                                        {        
                                                sum  = y->value * y->value;
                                                  y;
                                        }
                                        else
                                        {
                                                sum  = x->value * x->value;
                                                  x;
                                        }
                                }
                        }                            while(x->index != -1)
                        {
                                sum  = x->value * x->value;
                                  x;
                        }                            while(y->index != -1)
                        {
                                sum  = y->value * y->value;
                                  y;
                        }
                        
                        return exp(-param.gamma*sum);
                }
                case SIGMOID:
                        return tanh(param.gamma*dot(x,y) param.coef0);
                default:
                        return 0;        /* Unreachable */
        }
}    // 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:
        Solver() {};
        virtual ~Solver() {};            struct SolutionInfo {
                double obj;
                double rho;
                double upper_bound_p;
                double upper_bound_n;
                double r;        // for Solver_NU
        };            void Solve(int l, const QMatrix& Q, const double *b_, const schar *y_,
                   double *alpha_, double Cp, double Cn, double eps,
                   SolutionInfo* si, int shrinking);
protected:
        int active_size;
        schar *y;
        double *G;                // gradient of objective function
        enum { LOWER_BOUND, UPPER_BOUND, FREE };
        char *alpha_status;        // LOWER_BOUND, UPPER_BOUND, FREE
        double *alpha;
        const QMatrix *Q;
        const Qfloat *QD;
        double eps;
        double Cp,Cn;
        double *b;
        int *active_set;
        double *G_bar;                // gradient, if we treat free variables as 0
        int l;
        bool unshrinked;        // XXX            double get_C(int i)
        {
                return (y[i] > 0)? Cp : Cn;
        }
        void update_alpha_status(int i)
        {
                if(alpha[i] >= get_C(i))
                        alpha_status[i] = UPPER_BOUND;
                else if(alpha[i] <= 0)
                        alpha_status[i] = LOWER_BOUND;
                else alpha_status[i] = FREE;
        }
        bool is_upper_bound(int i) { return alpha_status[i] == UPPER_BOUND; }
        bool is_lower_bound(int i) { return alpha_status[i] == LOWER_BOUND; }
        bool is_free(int i) { return alpha_status[i] == FREE; }
        void swap_index(int i, int j);
        void reconstruct_gradient();
        virtual int select_working_set(int &i, int &j);
        virtual int max_violating_pair(int &i, int &j);
        virtual double calculate_rho();
        virtual void do_shrinking();
};    void Solver::swap_index(int i, int j)
{
        Q->swap_index(i,j);
        swap(y[i],y[j]);
        swap(G[i],G[j]);
        swap(alpha_status[i],alpha_status[j]);
        swap(alpha[i],alpha[j]);
        swap(b[i],b[j]);
        swap(active_set[i],active_set[j]);
        swap(G_bar[i],G_bar[j]);
}    void Solver::reconstruct_gradient()
{
        // reconstruct inactive elements of G from G_bar and free variables            if(active_size == l) return;            int i;
        for(i=active_size;iget_Q(i,l);
                        double alpha_i = alpha[i];
                        for(int j=active_size;jl = l;
        this->Q = &Q;
        QD=Q.get_QD();
        clone(b, b_,l);
        clone(y, y_,l);
        clone(alpha,alpha_,l);
        this->Cp = Cp;
        this->Cn = Cn;
        this->eps = eps;
        unshrinked = false;            // initialize alpha_status
        {
                alpha_status = new char[l];
                for(int i=0;i 0)
                        {
                                if(alpha[j] < 0)
                                {
                                        alpha[j] = 0;
                                        alpha[i] = diff;
                                }
                        }
                        else
                        {
                                if(alpha[i] < 0)
                                {
                                        alpha[i] = 0;
                                        alpha[j] = -diff;
                                }
                        }
                        if(diff > C_i - C_j)
                        {
                                if(alpha[i] > C_i)
                                {
                                        alpha[i] = C_i;
                                        alpha[j] = C_i - diff;
                                }
                        }
                        else
                        {
                                if(alpha[j] > C_j)
                                {
                                        alpha[j] = C_j;
                                        alpha[i] = C_j   diff;
                                }
                        }
                }
                else
                {
                        double quad_coef = Q_i[i] Q_j[j]-2*Q_i[j];
                        if (quad_coef <= 0)
                                quad_coef = TAU;
                        double delta = (G[i]-G[j])/quad_coef;
                        double sum = alpha[i]   alpha[j];
                        alpha[i] -= delta;
                        alpha[j]  = delta;                            if(sum > C_i)
                        {
                                if(alpha[i] > C_i)
                                {
                                        alpha[i] = C_i;
                                        alpha[j] = sum - C_i;
                                }
                        }
                        else
                        {
                                if(alpha[j] < 0)
                                {
                                        alpha[j] = 0;
                                        alpha[i] = sum;
                                }
                        }
                        if(sum > C_j)
                        {
                                if(alpha[j] > C_j)
                                {
                                        alpha[j] = C_j;
                                        alpha[i] = sum - C_j;
                                }
                        }
                        else
                        {
                                if(alpha[i] < 0)
                                {
                                        alpha[i] = 0;
                                        alpha[j] = sum;
                                }
                        }
                }                    // update G                    double delta_alpha_i = alpha[i] - old_alpha_i;
                double delta_alpha_j = alpha[j] - old_alpha_j;
                
                for(int k=0;krho = calculate_rho();            // calculate objective value
        {
                double v = 0;
                int i;
                for(i=0;iobj = v/2;
        }            // put back the solution
        {
                for(int i=0;iupper_bound_p = Cp;
        si->upper_bound_n = Cn;            info("\noptimization finished, #iter = %d\n",iter);            delete[] b;
        delete[] y;
        delete[] alpha;
        delete[] alpha_status;
        delete[] active_set;
        delete[] G;
        delete[] G_bar;
}    // return 1 if already optimal, return 0 otherwise
int Solver::select_working_set(int &out_i, int &out_j)
{
        // return i,j such that
        // i: maximizes -y_i * grad(f)_i, i in I_up(\alpha)
        // j: minimizes the decrease of obj value
        //    (if quadratic coefficeint <= 0, replace it with tau)
        //    -y_j*grad(f)_j < -y_i*grad(f)_i, j in I_low(\alpha)
        
        double Gmax = -INF;
        int Gmax_idx = -1;
        int Gmin_idx = -1;
        double obj_diff_min = INF;            for(int t=0;t= Gmax)
                                {
                                        Gmax = -G[t];
                                        Gmax_idx = t;
                                }
                }
                else
                {
                        if(!is_lower_bound(t))
                                if(G[t] >= Gmax)
                                {
                                        Gmax = G[t];
                                        Gmax_idx = t;
                                }
                }            int i = Gmax_idx;
        const Qfloat *Q_i = NULL;
        if(i != -1) // NULL Q_i not accessed: Gmax=-INF if i=-1
                Q_i = Q->get_Q(i,active_size);            for(int j=0;j= eps)
                                {
                                        double obj_diff; 
                                        double quad_coef=Q_i[i] QD[j]-2*y[i]*Q_i[j];
                                        if (quad_coef > 0)
                                                obj_diff = -(grad_diff*grad_diff)/quad_coef;
                                        else
                                                obj_diff = -(grad_diff*grad_diff)/TAU;                                            if (obj_diff <= obj_diff_min)
                                        {
                                                Gmin_idx=j;
                                                obj_diff_min = obj_diff;
                                        }
                                }
                        }
                }
                else
                {
                        if (!is_upper_bound(j))
                        {
                                double grad_diff= Gmax-G[j];
                                if (grad_diff >= eps)
                                {
                                        double obj_diff; 
                                        double quad_coef=Q_i[i] QD[j] 2*y[i]*Q_i[j];
                                        if (quad_coef > 0)
                                                obj_diff = -(grad_diff*grad_diff)/quad_coef;
                                        else
                                                obj_diff = -(grad_diff*grad_diff)/TAU;                                            if (obj_diff <= obj_diff_min)
                                        {
                                                Gmin_idx=j;
                                                obj_diff_min = obj_diff;
                                        }
                                }
                        }
                }
        }            if(Gmin_idx == -1)
                 return 1;            out_i = Gmax_idx;
        out_j = Gmin_idx;
        return 0;
}    // return 1 if already optimal, return 0 otherwise
int Solver::max_violating_pair(int &out_i, int &out_j)
{
        // return i,j: maximal violating pair            double Gmax1 = -INF;                // max { -y_i * grad(f)_i | i in I_up(\alpha) }
        int Gmax1_idx = -1;            double Gmax2 = -INF;                // max { y_i * grad(f)_i | i in I_low(\alpha) }
        int Gmax2_idx = -1;            for(int i=0;i= Gmax1)
                                {
                                        Gmax1 = -G[i];
                                        Gmax1_idx = i;
                                }
                        }
                        if(!is_lower_bound(i))        // d = -1
                        {
                                if(G[i] >= Gmax2)
                                {
                                        Gmax2 = G[i];
                                        Gmax2_idx = i;
                                }
                        }
                }
                else                // y = -1
                {
                        if(!is_upper_bound(i))        // d =  1
                        {
                                if(-G[i] >= Gmax2)
                                {
                                        Gmax2 = -G[i];
                                        Gmax2_idx = i;
                                }
                        }
                        if(!is_lower_bound(i))        // d = -1
                        {
                                if(G[i] >= Gmax1)
                                {
                                        Gmax1 = G[i];
                                        Gmax1_idx = i;
                                }
                        }
                }
        }            if(Gmax1 Gmax2 < eps)
                 return 1;            out_i = Gmax1_idx;
        out_j = Gmax2_idx;
        return 0;
}    void Solver::do_shrinking()
{
        int i,j,k;
        if(max_violating_pair(i,j)!=0) return;
        double Gm1 = -y[j]*G[j];
        double Gm2 = y[i]*G[i];            // shrink
        
        for(k=0;k= Gm1) continue;
                        }
                        else        if(-G[k] >= Gm2) continue;
                }
                else if(is_upper_bound(k))
                {
                        if(y[k]== 1)
                        {
                                if(G[k] >= Gm2) continue;
                        }
                        else        if(G[k] >= Gm1) continue;
                }
                else continue;                    --active_size;
                swap_index(k,active_size);
                --k;        // look at the newcomer
        }            // unshrink, check all variables again before final iterations            if(unshrinked || -(Gm1   Gm2) > eps*10) return;
        
        unshrinked = true;
        reconstruct_gradient();            for(k=l-1;k>=active_size;k--)
        {
                if(is_lower_bound(k))
                {
                        if(y[k]== 1)
                        {
                                if(-G[k] < Gm1) continue;
                        }
                        else        if(-G[k] < Gm2) continue;
                }
                else if(is_upper_bound(k))
                {
                        if(y[k]== 1)
                        {
                                if(G[k] < Gm2) continue;
                        }
                        else        if(G[k] < Gm1) continue;
                }
                else continue;                    swap_index(k,active_size);
                active_size  ;
                  k;        // look at the newcomer
        }
}    double Solver::calculate_rho()
{
        double r;
        int nr_free = 0;
        double ub = INF, lb = -INF, sum_free = 0;
        for(int i=0;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;
                }
        }            if(nr_free>0)
                r = sum_free/nr_free;
        else
                r = (ub lb)/2;            return r;
}    //
// Solver for nu-svm classification and regression
//
// additional constraint: e^T \alpha = constant
//
class Solver_NU : public Solver
{
public:
        Solver_NU() {}
        void Solve(int l, const QMatrix& Q, const double *b, const schar *y,
                   double *alpha, double Cp, double Cn, double eps,
                   SolutionInfo* si, int shrinking)
        {
                this->si = si;
                Solver::Solve(l,Q,b,y,alpha,Cp,Cn,eps,si,shrinking);
        }
private:
        SolutionInfo *si;
        int select_working_set(int &i, int &j);
        double calculate_rho();
        void do_shrinking();
};    // return 1 if already optimal, return 0 otherwise
int Solver_NU::select_working_set(int &out_i, int &out_j)
{
        // return i,j such that y_i = y_j and
        // i: maximizes -y_i * grad(f)_i, i in I_up(\alpha)
        // j: minimizes the decrease of obj value
        //    (if quadratic coefficeint <= 0, replace it with tau)
        //    -y_j*grad(f)_j < -y_i*grad(f)_i, j in I_low(\alpha)            double Gmaxp = -INF;
        int Gmaxp_idx = -1;            double Gmaxn = -INF;
        int Gmaxn_idx = -1;            int Gmin_idx = -1;
        double obj_diff_min = INF;            for(int t=0;t= Gmaxp)
                                {
                                        Gmaxp = -G[t];
                                        Gmaxp_idx = t;
                                }
                }
                else
                {
                        if(!is_lower_bound(t))
                                if(G[t] >= Gmaxn)
                                {
                                        Gmaxn = G[t];
                                        Gmaxn_idx = t;
                                }
                }            int ip = Gmaxp_idx;
        int in = Gmaxn_idx;
        const Qfloat *Q_ip = NULL;
        const Qfloat *Q_in = NULL;
        if(ip != -1) // NULL Q_ip not accessed: Gmaxp=-INF if ip=-1
                Q_ip = Q->get_Q(ip,active_size);
        if(in != -1)
                Q_in = Q->get_Q(in,active_size);            for(int j=0;j= eps)
                                {
                                        double obj_diff; 
                                        double quad_coef = Q_ip[ip] QD[j]-2*Q_ip[j];
                                        if (quad_coef > 0)
                                                obj_diff = -(grad_diff*grad_diff)/quad_coef;
                                        else
                                                obj_diff = -(grad_diff*grad_diff)/TAU;                                            if (obj_diff <= obj_diff_min)
                                        {
                                                Gmin_idx=j;
                                                obj_diff_min = obj_diff;
                                        }
                                }
                        }
                }
                else
                {
                        if (!is_upper_bound(j))
                        {
                                double grad_diff=Gmaxn-G[j];
                                if (grad_diff >= eps)
                                {
                                        double obj_diff; 
                                        double quad_coef = Q_in[in] QD[j]-2*Q_in[j];
                                        if (quad_coef > 0)
                                                obj_diff = -(grad_diff*grad_diff)/quad_coef;
                                        else
                                                obj_diff = -(grad_diff*grad_diff)/TAU;                                            if (obj_diff <= obj_diff_min)
                                        {
                                                Gmin_idx=j;
                                                obj_diff_min = obj_diff;
                                        }
                                }
                        }
                }
        }            if(Gmin_idx == -1)
                 return 1;            if (y[Gmin_idx] ==  1)
                out_i = Gmaxp_idx;
        else
                out_i = Gmaxn_idx;
        out_j = Gmin_idx;            return 0;
}    void Solver_NU::do_shrinking()
{
        double Gmax1 = -INF;        // max { -y_i * grad(f)_i | y_i =  1, i in I_up(\alpha) }
        double Gmax2 = -INF;        // max { y_i * grad(f)_i | y_i =  1, i in I_low(\alpha) }
        double Gmax3 = -INF;        // max { -y_i * grad(f)_i | y_i = -1, i in I_up(\alpha) }
        double Gmax4 = -INF;        // max { y_i * grad(f)_i | y_i = -1, i in I_low(\alpha) }            // find maximal violating pair first
        int k;
        for(k=0;k Gmax1) Gmax1 = -G[k];
                        }
                        else        if(-G[k] > Gmax3) Gmax3 = -G[k];
                }
                if(!is_lower_bound(k))
                {
                        if(y[k]== 1)
                        {        
                                if(G[k] > Gmax2) Gmax2 = G[k];
                        }
                        else        if(G[k] > Gmax4) Gmax4 = G[k];
                }
        }            // shrinking            double Gm1 = -Gmax2;
        double Gm2 = -Gmax1;
        double Gm3 = -Gmax4;
        double Gm4 = -Gmax3;            for(k=0;k= Gm1) continue;
                        }
                        else        if(-G[k] >= Gm3) continue;
                }
                else if(is_upper_bound(k))
                {
                        if(y[k]== 1)
                        {
                                if(G[k] >= Gm2) continue;
                        }
                        else        if(G[k] >= Gm4) continue;
                }
                else continue;                    --active_size;
                swap_index(k,active_size);
                --k;        // look at the newcomer
        }            // unshrink, check all variables again before final iterations            if(unshrinked || max(-(Gm1 Gm2),-(Gm3 Gm4)) > eps*10) return;
        
        unshrinked = true;
        reconstruct_gradient();            for(k=l-1;k>=active_size;k--)
        {
                if(is_lower_bound(k))
                {
                        if(y[k]== 1)
                        {
                                if(-G[k] < Gm1) continue;
                        }
                        else        if(-G[k] < Gm3) continue;
                }
                else if(is_upper_bound(k))
                {
                        if(y[k]== 1)
                        {
                                if(G[k] < Gm2) continue;
                        }
                        else        if(G[k] < Gm4) continue;
                }
                else continue;                    swap_index(k,active_size);
                active_size  ;
                  k;        // look at the newcomer
        }
}    double Solver_NU::calculate_rho()
{
        int nr_free1 = 0,nr_free2 = 0;
        double ub1 = INF, ub2 = INF;
        double lb1 = -INF, lb2 = -INF;
        double sum_free1 = 0, sum_free2 = 0;            for(int i=0;i 0)
                r1 = sum_free1/nr_free1;
        else
                r1 = (ub1 lb1)/2;
        
        if(nr_free2 > 0)
                r2 = sum_free2/nr_free2;
        else
                r2 = (ub2 lb2)/2;
        
        si->r = (r1 r2)/2;
        return (r1-r2)/2;
}    //
// Q matrices for various formulations
//
class SVC_Q: public Kernel
{ 
public:
        SVC_Q(const svm_problem& prob, const svm_parameter& param, const schar *y_)
        :Kernel(prob.l, prob.x, param)
        {
                clone(y,y_,prob.l);
                cache = new Cache(prob.l,(int)(param.cache_size*(1<<20)));
                QD = new Qfloat[prob.l];
                for(int i=0;i*kernel_function)(i,i);
        }
        
        Qfloat *get_Q(int i, int len) const
        {
                Qfloat *data;
                int start;
                if((start = cache->get_data(i,&data,len)) < len)
                {
                        for(int j=start;j*kernel_function)(i,j));
                }
                return data;
        }            Qfloat *get_QD() const
        {
                return QD;
        }            void swap_index(int i, int j) const
        {
                cache->swap_index(i,j);
                Kernel::swap_index(i,j);
                swap(y[i],y[j]);
                swap(QD[i],QD[j]);
        }            ~SVC_Q()
        {
                delete[] y;
                delete cache;
                delete[] QD;
        }
private:
        schar *y;
        Cache *cache;
        Qfloat *QD;
};    class ONE_CLASS_Q: public Kernel
{
public:
        ONE_CLASS_Q(const svm_problem& prob, const svm_parameter& param)
        :Kernel(prob.l, prob.x, param)
        {
                cache = new Cache(prob.l,(int)(param.cache_size*(1<<20)));
                QD = new Qfloat[prob.l];
                for(int i=0;i*kernel_function)(i,i);
        }
        
        Qfloat *get_Q(int i, int len) const
        {
                Qfloat *data;
                int start;
                if((start = cache->get_data(i,&data,len)) < len)
                {
                        for(int j=start;j*kernel_function)(i,j);
                }
                return data;
        }            Qfloat *get_QD() const
        {
                return QD;
        }            void swap_index(int i, int j) const
        {
                cache->swap_index(i,j);
                Kernel::swap_index(i,j);
                swap(QD[i],QD[j]);
        }            ~ONE_CLASS_Q()
        {
                delete cache;
                delete[] QD;
        }
private:
        Cache *cache;
        Qfloat *QD;
};    class SVR_Q: public Kernel
{ 
public:
        SVR_Q(const svm_problem& prob, const svm_parameter& param)
        :Kernel(prob.l, prob.x, param)
        {
                l = prob.l;
                cache = new Cache(l,(int)(param.cache_size*(1<<20)));
                QD = new Qfloat[2*l];
                sign = new schar[2*l];
                index = new int[2*l];
                for(int k=0;k*kernel_function)(k,k);
                        QD[k l]=QD[k];
                }
                buffer[0] = new Qfloat[2*l];
                buffer[1] = new Qfloat[2*l];
                next_buffer = 0;
        }            void swap_index(int i, int j) const
        {
                swap(sign[i],sign[j]);
                swap(index[i],index[j]);
                swap(QD[i],QD[j]);
        }
        
        Qfloat *get_Q(int i, int len) const
        {
                Qfloat *data;
                int real_i = index[i];
                if(cache->get_data(real_i,&data,l) < l)
                {
                        for(int j=0;j*kernel_function)(real_i,j);
                }                    // reorder and copy
                Qfloat *buf = buffer[next_buffer];
                next_buffer = 1 - next_buffer;
                schar si = sign[i];
                for(int j=0;jl;
        double *minus_ones = new double[l];
        schar *y = new schar[l];            int i;            for(i=0;iy[i] > 0) y[i] =  1; else y[i]=-1;
        }            Solver s;
        s.Solve(l, SVC_Q(*prob,*param,y), minus_ones, y,
                alpha, Cp, Cn, param->eps, si, param->shrinking);            double sum_alpha=0;
        for(i=0;il));            for(i=0;il;
        double nu = param->nu;            schar *y = new schar[l];            for(i=0;iy[i]>0)
                        y[i] =  1;
                else
                        y[i] = -1;            double sum_pos = nu*l/2;
        double sum_neg = nu*l/2;            for(i=0;ieps, si,  param->shrinking);
        double r = si->r;            info("C = %f\n",1/r);            for(i=0;irho /= r;
        si->obj /= (r*r);
        si->upper_bound_p = 1/r;
        si->upper_bound_n = 1/r;            delete[] y;
        delete[] zeros;
}    static void solve_one_class(
        const svm_problem *prob, const svm_parameter *param,
        double *alpha, Solver::SolutionInfo* si)
{
        int l = prob->l;
        double *zeros = new double[l];
        schar *ones = new schar[l];
        int i;            int n = (int)(param->nu*prob->l);        // # of alpha's at upper bound            for(i=0;inu * prob->l - n;
        for(i=n 1;ieps, si, param->shrinking);            delete[] zeros;
        delete[] ones;
}    static void solve_epsilon_svr(
        const svm_problem *prob, const svm_parameter *param,
        double *alpha, Solver::SolutionInfo* si)
{
        int l = prob->l;
        double *alpha2 = new double[2*l];
        double *linear_term = new double[2*l];
        schar *y = new schar[2*l];
        int i;            for(i=0;ip - prob->y[i];
                y[i] = 1;                    alpha2[i l] = 0;
                linear_term[i l] = param->p   prob->y[i];
                y[i l] = -1;
        }            Solver s;
        s.Solve(2*l, SVR_Q(*prob,*param), linear_term, y,
                alpha2, param->C, param->C, param->eps, si, param->shrinking);            double sum_alpha = 0;
        for(i=0;iC*l));            delete[] alpha2;
        delete[] linear_term;
        delete[] y;
}    static void solve_nu_svr(
        const svm_problem *prob, const svm_parameter *param,
        double *alpha, Solver::SolutionInfo* si)
{
        int l = prob->l;
        double C = param->C;
        double *alpha2 = new double[2*l];
        double *linear_term = new double[2*l];
        schar *y = new schar[2*l];
        int i;            double sum = C * param->nu * l / 2;
        for(i=0;iy[i];
                y[i] = 1;                    linear_term[i l] = prob->y[i];
                y[i l] = -1;
        }            Solver_NU s;
        s.Solve(2*l, SVR_Q(*prob,*param), linear_term, y,
                alpha2, C, C, param->eps, si, param->shrinking);            info("epsilon = %f\n",-si->r);            for(i=0;il);
        Solver::SolutionInfo si;
        switch(param->svm_type)
        {
                case C_SVC:
                        solve_c_svc(prob,param,alpha,&si,Cp,Cn);
                        break;
                case NU_SVC:
                        solve_nu_svc(prob,param,alpha,&si);
                        break;
                case ONE_CLASS:
                        solve_one_class(prob,param,alpha,&si);
                        break;
                case EPSILON_SVR:
                        solve_epsilon_svr(prob,param,alpha,&si);
                        break;
                case NU_SVR:
                        solve_nu_svr(prob,param,alpha,&si);
                        break;
        }            info("obj = %f, rho = %f\n",si.ob