svm.cpp

/*
Copyright (c) 2000-2007 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.
*/

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <ctype.h>
#include <float.h>
#include <string.h>
#include <stdarg.h>
#include "svm.h"
typedef float Qfloat;
typedef signed char schar;
#ifndef min
template <class T> inline T min(T x,T y) { return (x<y)?x:y; }
#endif
#ifndef max
template <class T> inline T max(T x,T y) { return (x>y)?x:y; }
#endif
template <class T> inline void swap(T& x, T& y) { T t=x; x=y; y=t; }
template <class S, class T> inline void clone(T*& dst, S* src, int n)
{
        dst = new T[n];
        memcpy((void *)dst,(void *)src,sizeof(T)*n);
}
inline double powi(double base, int times)
{
        double tmp = base, ret = 1.0;

        for(int t=times; t>0; t/=2)
        {
                if(t%2==1) ret*=tmp;
                tmp = tmp * tmp;
        }
        return ret;
}
#define INF HUGE_VAL
#define TAU 1e-12
#define Malloc(type,n) (type *)malloc((n)*sizeof(type))
#if 1
void info(const 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,long 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;
        long 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_,long 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 * (long int) 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;
        virtual ~QMatrix() {}
};

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 int 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 powi(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);
        }
        double kernel_precomputed(int i, int j) const
        {
                return x[i][(int)(x[j][0].value)].value;
        }
};

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;
                case PRECOMPUTED:
                        kernel_function = &Kernel::kernel_precomputed;
                        break;
        }

        clone(x,x_,l);

        if(kernel_type == RBF)
        {
                x_square = new double[l];
                for(int i=0;i<l;i++)
                        x_square[i] = dot(x[i],x[i]);
        }
        else
                x_square = 0;
}

Kernel::~Kernel()
{
        delete[] x;
        delete[] x_square;
}

double Kernel::dot(const svm_node *px, const svm_node *py)
{
        double sum = 0;
        while(px->index != -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 powi(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);
                case PRECOMPUTED:  //x: test (validation), y: SV
                        return x[(int)(y->value)].value;
                default:
                        return 0;  // Unreachable
        }
}

// An SMO algorithm in Fan et al., JMLR 6(2005), p. 1889--1918
// Solves:
//
//        min 0.5(\alpha^T Q \alpha) + p^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, p, y, Cp, Cn, and an initial feasible point \alpha
//        l is the size of vectors and matrices
//        eps is the stopping tolerance
//
// 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 *p_, 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 *p;
        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 double calculate_rho();
        virtual void do_shrinking();
private:
        bool be_shrunken(int i, double Gmax1, double Gmax2);
};

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(p[i],p[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;i<l;i++)
                G[i] = G_bar[i] + p[i];

        for(i=0;i<active_size;i++)
                if(is_free(i))
                {
                        const Qfloat *Q_i = Q->get_Q(i,l);
                        double alpha_i = alpha[i];
                        for(int j=active_size;j<l;j++)
                                G[j] += alpha_i * Q_i[j];
                }
}

void Solver::Solve(int l, const QMatrix& Q, const double *p_, const schar *y_,
                   double *alpha_, double Cp, double Cn, double eps,
                   SolutionInfo* si, int shrinking)
{
        this->l = l;
        this->Q = &Q;
        QD=Q.get_QD();
        clone(p, p_,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<l;i++)
                        update_alpha_status(i);
        }

        // initialize active set (for shrinking)
        {
                active_set = new int[l];
                for(int i=0;i<l;i++)
                        active_set[i] = i;
                active_size = l;
        }

        // initialize gradient
        {
                G = new double[l];
                G_bar = new double[l];
                int i;
                for(i=0;i<l;i++)
                {
                        G[i] = p[i];
                        G_bar[i] = 0;
                }
                for(i=0;i<l;i++)
                        if(!is_lower_bound(i))
                        {
                                const Qfloat *Q_i = Q.get_Q(i,l);
                                double alpha_i = alpha[i];
                                int j;
                                for(j=0;j<l;j++)
                                        G[j] += alpha_i*Q_i[j];
                                if(is_upper_bound(i))
                                        for(j=0;j<l;j++)
                                                G_bar[j] += get_C(i) * Q_i[j];
                        }
        }

        // optimization step

        int iter = 0;
        int counter = min(l,1000)+1;

        while(1)
        {
                // show progress and do shrinking

                if(--counter == 0)
                {
                        counter = min(l,1000);
                        if(shrinking) do_shrinking();
                        info("."); info_flush();
                }

                int i,j;
                if(select_working_set(i,j)!=0)
                {
                        // reconstruct the whole gradient
                        reconstruct_gradient();
                        // reset active set size and check
                        active_size = l;
                        info("*"); info_flush();
                        if(select_working_set(i,j)!=0)
                                break;
                        else
                                counter = 1;        // do shrinking next iteration
                }

                ++iter;

                // update alpha[i] and alpha[j], handle bounds carefully

                const Qfloat *Q_i = Q.get_Q(i,active_size);
                const Qfloat *Q_j = Q.get_Q(j,active_size);

                double C_i = get_C(i);
                double C_j = get_C(j);

                double old_alpha_i = alpha[i];
                double old_alpha_j = alpha[j];

                if(y[i]!=y[j])
                {
                        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 diff = alpha[i] - alpha[j];
                        alpha[i] += delta;
                        alpha[j] += delta;

                        if(diff > 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;k<active_size;k++)
                {
                        G[k] += Q_i[k]*delta_alpha_i + Q_j[k]*delta_alpha_j;
                }

                // update alpha_status and G_bar

                {
                        bool ui = is_upper_bound(i);
                        bool uj = is_upper_bound(j);
                        update_alpha_status(i);
                        update_alpha_status(j);
                        int k;
                        if(ui != is_upper_bound(i))
                        {
                                Q_i = Q.get_Q(i,l);
                                if(ui)
                                        for(k=0;k<l;k++)
                                                G_bar[k] -= C_i * Q_i[k];
                                else
                                        for(k=0;k<l;k++)
                                                G_bar[k] += C_i * Q_i[k];
                        }

                        if(uj != is_upper_bound(j))
                        {
                                Q_j = Q.get_Q(j,l);
                                if(uj)
                                        for(k=0;k<l;k++)
                                                G_bar[k] -= C_j * Q_j[k];
                                else
                                        for(k=0;k<l;k++)
                                                G_bar[k] += C_j * Q_j[k];
                        }
                }
        }

        // calculate rho

        si->rho = calculate_rho();

        // calculate objective value
        {
                double v = 0;
                int i;
                for(i=0;i<l;i++)
                        v += alpha[i] * (G[i] + p[i]);

                si->obj = v/2;
        }

        // put back the solution
        {
                for(int i=0;i<l;i++)
                        alpha_[active_set[i]] = alpha[i];
        }

        // juggle everything back
        /*{
                for(int i=0;i<l;i++)
                        while(active_set[i] != i)
                                swap_index(i,active_set[i]);
                                // or Q.swap_index(i,active_set[i]);
        }*/

        si->upper_bound_p = Cp;
        si->upper_bound_n = Cn;

        info("\noptimization finished, #iter = %d\n",iter);

        delete[] p;
        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;
        double Gmax2 = -INF;
        int Gmax_idx = -1;
        int Gmin_idx = -1;
        double obj_diff_min = INF;

        for(int t=0;t<active_size;t++)
                if(y[t]==+1)
                {
                        if(!is_upper_bound(t))
                                if(-G[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<active_size;j++)
        {
                if(y[j]==+1)
                {
                        if (!is_lower_bound(j))
                        {
                                double grad_diff=Gmax+G[j];
                                if (G[j] >= Gmax2)
                                        Gmax2 = G[j];
                                if (grad_diff > 0)
                                {
                                        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 (-G[j] >= Gmax2)
                                        Gmax2 = -G[j];
                                if (grad_diff > 0)
                                {
                                        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(Gmax+Gmax2 < eps)
                return 1;

        out_i = Gmax_idx;
        out_j = Gmin_idx;
        return 0;
}

bool Solver::be_shrunken(int i, double Gmax1, double Gmax2)
{
        if(is_upper_bound(i))
        {
                if(y[i]==+1)
                        return(-G[i] > Gmax1);
                else
                        return(-G[i] > Gmax2);
        }
        else if(is_lower_bound(i))
        {
                if(y[i]==+1)
                        return(G[i] > Gmax2);
                else
                        return(G[i] > Gmax1);
        }
        else
                return(false);
}

void Solver::do_shrinking()
{
        int i;
        double Gmax1 = -INF;                // max { -y_i * grad(f)_i | i in I_up(\alpha) }
        double Gmax2 = -INF;                // max { y_i * grad(f)_i | i in I_low(\alpha) }

        // find maximal violating pair first
        for(i=0;i<active_size;i++)
        {
                if(y[i]==+1)
                {
                        if(!is_upper_bound(i))
                        {
                                if(-G[i] >= Gmax1)
                                        Gmax1 = -G[i];
                        }
                        if(!is_lower_bound(i))
                        {
                                if(G[i] >= Gmax2)
                                        Gmax2 = G[i];
                        }
                }
                else
                {
                        if(!is_upper_bound(i))
                        {
                                if(-G[i] >= Gmax2)
                                        Gmax2 = -G[i];
                        }
                        if(!is_lower_bound(i))
                        {
                                if(G[i] >= Gmax1)
                                        Gmax1 = G[i];
                        }
                }
        }

        // shrink

        for(i=0;i<active_size;i++)
                if (be_shrunken(i, Gmax1, Gmax2))
                {
                        active_size--;
                        while (active_size > i)
                        {
                                if (!be_shrunken(active_size, Gmax1, Gmax2))
                                {
                                        swap_index(i,active_size);
                                        break;
                                }
                                active_size--;
                        }
                }

        // unshrink, check all variables again before final iterations

        if(unshrinked || Gmax1 + Gmax2 > eps*10) return;

        unshrinked = true;
        reconstruct_gradient();

        for(i=l-1;i>=active_size;i--)
                if (!be_shrunken(i, Gmax1, Gmax2))
                {
                        while (active_size < i)
                        {
                                if (be_shrunken(active_size, Gmax1, Gmax2))
                                {
                                        swap_index(i,active_size);
                                        break;
                                }
                                active_size++;
                        }
                        active_size++;
                }
}

double Solver::calculate_rho()
{
        double r;
        int nr_free = 0;
        double ub = INF, lb = -INF, sum_free = 0;
        for(int i=0;i<active_size;i++)
        {
                double yG = y[i]*G[i];

                if(is_upper_bound(i))
                {
                        if(y[i]==-1)
                                ub = min(ub,yG);
                        else
                                lb = max(lb,yG);
                }
                else if(is_lower_bound(i))
                {
                        if(y[i]==+1)
                                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 *p, const schar *y,
                   double *alpha, double Cp, double Cn, double eps,
                   SolutionInfo* si, int shrinking)
        {
                this->si = si;
                Solver::Solve(l,Q,p,y,alpha,Cp,Cn,eps,si,shrinking);
        }
private:
        SolutionInfo *si;
        int select_working_set(int &i, int &j);
        double calculate_rho();
        bool be_shrunken(int i, double Gmax1, double Gmax2, double Gmax3, double Gmax4);
        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;
        double Gmaxp2 = -INF;
        int Gmaxp_idx = -1;

        double Gmaxn = -INF;
        double Gmaxn2 = -INF;
        int Gmaxn_idx = -1;

        int Gmin_idx = -1;
        double obj_diff_min = INF;

        for(int t=0;t<active_size;t++)
                if(y[t]==+1)
                {
                        if(!is_upper_bound(t))
                                if(-G[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<active_size;j++)
        {
                if(y[j]==+1)
                {
                        if (!is_lower_bound(j))
                        {
                                double grad_diff=Gmaxp+G[j];
                                if (G[j] >= Gmaxp2)
                                        Gmaxp2 = G[j];
                                if (grad_diff > 0)
                                {
                                        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 (-G[j] >= Gmaxn2)
                                        Gmaxn2 = -G[j];
                                if (grad_diff > 0)
                                {
                                        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(max(Gmaxp+Gmaxp2,Gmaxn+Gmaxn2) < eps)
                 return 1;

        if (y[Gmin_idx] == +1)
                out_i = Gmaxp_idx;
        else
                out_i = Gmaxn_idx;
        out_j = Gmin_idx;

        return 0;
}

bool Solver_NU::be_shrunken(int i, double Gmax1, double Gmax2, double Gmax3, double Gmax4)
{
        if(is_upper_bound(i))
        {
                if(y[i]==+1)
                        return(-G[i] > Gmax1);
                else
                        return(-G[i] > Gmax4);
        }
        else if(is_lower_bound(i))
        {
                if(y[i]==+1)
                        return(G[i] > Gmax2);
                else
                        return(G[i] > Gmax3);
        }
        else
                return(false);
}

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 i;
        for(i=0;i<active_size;i++)
        {
                if(!is_upper_bound(i))
                {
                        if(y[i]==+1)
                        {
                                if(-G[i] > Gmax1) Gmax1 = -G[i];
                        }
                        else        if(-G[i] > Gmax4) Gmax4 = -G[i];
                }
                if(!is_lower_bound(i))
                {
                        if(y[i]==+1)
                        {
                                if(G[i] > Gmax2) Gmax2 = G[i];
                        }
                        else        if(G[i] > Gmax3) Gmax3 = G[i];
                }
        }

        // shrinking

        for(i=0;i<active_size;i++)
                if (be_shrunken(i, Gmax1, Gmax2, Gmax3, Gmax4))
                {
                        active_size--;
                        while (active_size > i)
                        {
                                if (!be_shrunken(active_size, Gmax1, Gmax2, Gmax3, Gmax4))
                                {
                                        swap_index(i,active_size);
                                        break;
                                }
                                active_size--;
                        }
                }

        // unshrink, check all variables again before final iterations

        if(unshrinked || max(Gmax1+Gmax2,Gmax3+Gmax4) > eps*10) return;

        unshrinked = true;
        reconstruct_gradient();

        for(i=l-1;i>=active_size;i--)
                if (!be_shrunken(i, Gmax1, Gmax2, Gmax3, Gmax4))
                {
                        while (active_size < i)
                        {
                                if (be_shrunken(active_size, Gmax1, Gmax2, Gmax3, Gmax4))
                                {
                                        swap_index(i,active_size);
                                        break;
                                }
                                active_size++;
                        }
                        active_size++;
                }
}

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<active_size;i++)
        {
                if(y[i]==+1)
                {
                        if(is_upper_bound(i))
                                lb1 = max(lb1,G[i]);
                        else if(is_lower_bound(i))
                                ub1 = min(ub1,G[i]);
                        else
                        {
                                ++nr_free1;
                                sum_free1 += G[i];
                        }
                }
                else
                {
                        if(is_upper_bound(i))
                                lb2 = max(lb2,G[i]);
                        else if(is_lower_bound(i))
                                ub2 = min(ub2,G[i]);
                        else
                        {
                                ++nr_free2;
                                sum_free2 += G[i];
                        }
                }
        }

        double r1,r2;
        if(nr_free1 > 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,(long int)(param.cache_size*(1<<20)));
                QD = new Qfloat[prob.l];
                for(int i=0;i<prob.l;i++)
                        QD[i]= (Qfloat)(this->*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<len;j++)
                                data[j] = (Qfloat)(y[i]*y[j]*(this->*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,(long int)(param.cache_size*(1<<20)));
                QD = new Qfloat[prob.l];
                for(int i=0;i<prob.l;i++)
                        QD[i]= (Qfloat)(this->*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<len;j++)
                                data[j] = (Qfloat)(this->*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,(long 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<l;k++)
                {
                        sign[k] = 1;
                        sign[k+l] = -1;
                        index[k] = k;
                        index[k+l] = k;
                        QD[k]= (Qfloat)(this->*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<l;j++)
                                data[j] = (Qfloat)(this->*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;j<len;j++)
                        buf[j] = si * sign[j] * data[index[j]];
                return buf;
        }

        Qfloat *get_QD() const
        {
                return QD;
        }

        ~SVR_Q()
        {
                delete cache;
                delete[] sign;
                delete[] index;
                delete[] buffer[0];
                delete[] buffer[1];
                delete[] QD;
        }
private:
        int l;
        Cache *cache;
        schar *sign;
        int *index;
        mutable int next_buffer;
        Qfloat *buffer[2];
        Qfloat *QD;
};

//
// construct and solve various formulations
//
static void solve_c_svc(
        const svm_problem *prob, const svm_parameter* param,
        double *alpha, Solver::SolutionInfo* si, double Cp, double Cn)
{
        int l = prob->l;
        double *minus_ones = new double[l];
        schar *y = new schar[l];

        int i;

        for(i=0;i<l;i++)
        {
                alpha[i] = 0;
                minus_ones[i] = -1;
                if(prob->y[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;i<l;i++)
                sum_alpha += alpha[i];

        if (Cp==Cn)
                info("nu = %f\n", sum_alpha/(Cp*prob->l));

        for(i=0;i<l;i++)
                alpha[i] *= y[i];

        delete[] minus_ones;
        delete[] y;
}

static void solve_nu_svc(
        const svm_problem *prob, const svm_parameter *param,
        double *alpha, Solver::SolutionInfo* si)
{
        int i;
        int l = prob->l;
        double nu = param->nu;

        schar *y = new schar[l];

        for(i=0;i<l;i++)
                if(prob->y[i]>0)
                        y[i] = +1;
                else
                        y[i] = -1;

        double sum_pos = nu*l/2;
        double sum_neg = nu*l/2;

        for(i=0;i<l;i++)
                if(y[i] == +1)
                {
                        alpha[i] = min(1.0,sum_pos);
                        sum_pos -= alpha[i];
                }
                else
                {
                        alpha[i] = min(1.0,sum_neg);
                        sum_neg -= alpha[i];
                }

        double *zeros = new double[l];

        for(i=0;i<l;i++)
                zeros[i] = 0;

        Solver_NU s;
        s.Solve(l, SVC_Q(*prob,*param,y), zeros, y,
                alpha, 1.0, 1.0, param->eps, si,  param->shrinking);
        double r = si->r;

        info("C = %f\n",1/r);

        for(i=0;i<l;i++)
                alpha[i] *= y[i]/r;

        si->rho /= 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;i<n;i++)
                alpha[i] = 1;
        if(n<prob->l)
                alpha[n] = param->nu * prob->l - n;
        for(i=n+1;i<l;i++)
                alpha[i] = 0;

        for(i=0;i<l;i++)
        {
                zeros[i] = 0;
                ones[i] = 1;
        }

        Solver s;
        s.Solve(l, ONE_CLASS_Q(*prob,*param), zeros, ones,
                alpha, 1.0, 1.0, param->eps, 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;i<l;i++)
        {
                alpha2[i] = 0;
                linear_term[i] = param->p - 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;i<l;i++)
        {
                alpha[i] = alpha2[i] - alpha2[i+l];
                sum_alpha += fabs(alpha[i]);
        }
        info("nu = %f\n",sum_alpha/(param->C*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;i<l;i++)
        {
                alpha2[i] = alpha2[i+l] = min(sum,C);
                sum -= alpha2[i];

                linear_term[i] = - prob->y[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;i<l;i++)
                alpha[i] = alpha2[i] - alpha2[i+l];

        delete[] alpha2;
        delete[] linear_term;
        delete[] y;
}

//
// decision_function
//
struct decision_function
{
        double *alpha;
        double rho;
};

decision_function svm_train_one(
        const svm_problem *prob, const svm_parameter *param,
        double Cp, double Cn)
{
        double *alpha = Malloc(double,prob->l);
        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.obj,si.rho);

        // output SVs

        int nSV = 0;
        int nBSV = 0;
        for(int i=0;i<prob->l;i++)
        {
                if(fabs(alpha[i]) > 0)
                {
                        ++nSV;
                        if(prob->y[i] > 0)
                        {
                                if(fabs(alpha[i]) >= si.upper_bound_p)
                                        ++nBSV;
                        }
                        else
                        {
                                if(fabs(alpha[i]) >= si.upper_bound_n)
                                        ++nBSV;
                        }
                }
        }

        info("nSV = %d, nBSV = %d\n",nSV,nBSV);

        decision_function f;
        f.alpha = alpha;
        f.rho = si.rho;
        return f;
}

//
// svm_model
//
struct svm_model
{
        svm_parameter param;        // parameter
        int nr_class;                // number of classes, = 2 in regression/one class svm
        int l;                        // total #SV
        svm_node **SV;                // SVs (SV[l])
        double **sv_coef;        // coefficients for SVs in decision functions (sv_coef[k-1][l])
        double *rho;                // constants in decision functions (rho[k*(k-1)/2])
        double *probA;          // pariwise probability information
        double *probB;

        // for classification only

        int *label;                // label of each class (label[k])
        int *nSV;                // number of SVs for each class (nSV[k])
                                // nSV[0] + nSV[1] + ... + nSV[k-1] = l
        // XXX
        int free_sv;                // 1 if svm_model is created by svm_load_model
                                // 0 if svm_model is created by svm_train
};

// Platt's binary SVM Probablistic Output: an improvement from Lin et al.
void sigmoid_train(
        int l, const double *dec_values, const double *labels,
        double& A, double& B)
{
        double prior1=0, prior0 = 0;
        int i;

        for (i=0;i<l;i++)
                if (labels[i] > 0) prior1+=1;
                else prior0+=1;

        int max_iter=100;         // Maximal number of iterations
        double min_step=1e-10;        // Minimal step taken in line search
        double sigma=1e-12;        // For numerically strict PD of Hessian
        double eps=1e-5;
        double hiTarget=(prior1+1.0)/(prior1+2.0);
        double loTarget=1/(prior0+2.0);
        double *t=Malloc(double,l);
        double fApB,p,q,h11,h22,h21,g1,g2,det,dA,dB,gd,stepsize;
        double newA,newB,newf,d1,d2;
        int iter;

        // Initial Point and Initial Fun Value
        A=0.0; B=log((prior0+1.0)/(prior1+1.0));
        double fval = 0.0;

        for (i=0;i<l;i++)
        {
                if (labels[i]>0) t[i]=hiTarget;
                else t[i]=loTarget;
                fApB = dec_values[i]*A+B;
                if (fApB>=0)
                        fval += t[i]*fApB + log(1+exp(-fApB));
                else
                        fval += (t[i] - 1)*fApB +log(1+exp(fApB));
        }
        for (iter=0;iter<max_iter;iter++)
        {
                // Update Gradient and Hessian (use H' = H + sigma I)
                h11=sigma; // numerically ensures strict PD
                h22=sigma;
                h21=0.0;g1=0.0;g2=0.0;
                for (i=0;i<l;i++)
                {
                        fApB = dec_values[i]*A+B;
                        if (fApB >= 0)
                        {
                                p=exp(-fApB)/(1.0+exp(-fApB));
                                q=1.0/(1.0+exp(-fApB));
                        }
                        else
                        {
                                p=1.0/(1.0+exp(fApB));
                                q=exp(fApB)/(1.0+exp(fApB));
                        }
                        d2=p*q;
                        h11+=dec_values[i]*dec_values[i]*d2;
                        h22+=d2;
                        h21+=dec_values[i]*d2;
                        d1=t[i]-p;
                        g1+=dec_values[i]*d1;
                        g2+=d1;
                }

                // Stopping Criteria
                if (fabs(g1)<eps && fabs(g2)<eps)
                        break;

                // Finding Newton direction: -inv(H') * g
                det=h11*h22-h21*h21;
                dA=-(h22*g1 - h21 * g2) / det;
                dB=-(-h21*g1+ h11 * g2) / det;
                gd=g1*dA+g2*dB;


                stepsize = 1;                 // Line Search
                while (stepsize >= min_step)
                {
                        newA = A + stepsize * dA;
                        newB = B + stepsize * dB;

                        // New function value
                        newf = 0.0;
                        for (i=0;i<l;i++)
                        {
                                fApB = dec_values[i]*newA+newB;
                                if (fApB >= 0)
                                        newf += t[i]*fApB + log(1+exp(-fApB));
                                else
                                        newf += (t[i] - 1)*fApB +log(1+exp(fApB));
                        }
                        // Check sufficient decrease
                        if (newf<fval+0.0001*stepsize*gd)
                        {
                                A=newA;B=newB;fval=newf;
                                break;
                        }
                        else
                                stepsize = stepsize / 2.0;
                }

                if (stepsize < min_step)
                {
                        info("Line search fails in two-class probability estimates\n");
                        break;
                }
        }

        if (iter>=max_iter)
                info("Reaching maximal iterations in two-class probability estimates\n");
        free(t);
}

double sigmoid_predict(double decision_value, double A, double B)
{
        double fApB = decision_value*A+B;
        if (fApB >= 0)
                return exp(-fApB)/(1.0+exp(-fApB));
        else
                return 1.0/(1+exp(fApB)) ;
}

// Method 2 from the multiclass_prob paper by Wu, Lin, and Weng
void multiclass_probability(int k, double **r, double *p)
{
        int t,j;
        int iter = 0, max_iter=max(100,k);
        double **Q=Malloc(double *,k);
        double *Qp=Malloc(double,k);
        double pQp, eps=0.005/k;

        for (t=0;t<k;t++)
        {
                p[t]=1.0/k;  // Valid if k = 1
                Q[t]=Malloc(double,k);
                Q[t][t]=0;
                for (j=0;j<t;j++)
                {
                        Q[t][t]+=r[j][t]*r[j][t];
                        Q[t][j]=Q[j][t];
                }
                for (j=t+1;j<k;j++)
                {
                        Q[t][t]+=r[j][t]*r[j][t];
                        Q[t][j]=-r[j][t]*r[t][j];
                }
        }
        for (iter=0;iter<max_iter;iter++)
        {
                // stopping condition, recalculate QP,pQP for numerical accuracy
                pQp=0;
                for (t=0;t<k;t++)
                {
                        Qp[t]=0;
                        for (j=0;j<k;j++)
                                Qp[t]+=Q[t][j]*p[j];
                        pQp+=p[t]*Qp[t];
                }
                double max_error=0;
                for (t=0;t<k;t++)
                {
                        double error=fabs(Qp[t]-pQp);
                        if (error>max_error)
                                max_error=error;
                }
                if (max_error<eps) break;

                for (t=0;t<k;t++)
                {
                        double diff=(-Qp[t]+pQp)/Q[t][t];
                        p[t]+=diff;
                        pQp=(pQp+diff*(diff*Q[t][t]+2*Qp[t]))/(1+diff)/(1+diff);
                        for (j=0;j<k;j++)
                        {
                                Qp[j]=(Qp[j]+diff*Q[t][j])/(1+diff);
                                p[j]/=(1+diff);
                        }
                }
        }
        if (iter>=max_iter)
                info("Exceeds max_iter in multiclass_prob\n");
        for(t=0;t<k;t++) free(Q[t]);
        free(Q);
        free(Qp);
}

// Cross-validation decision values for probability estimates
void svm_binary_svc_probability(
        const svm_problem *prob, const svm_parameter *param,
        double Cp, double Cn, double& probA, double& probB)
{
        int i;
        int nr_fold = 5;
        int *perm = Malloc(int,prob->l);
        double *dec_values = Malloc(double,prob->l);

        // random shuffle
        for(i=0;i<prob->l;i++) perm[i]=i;
        for(i=0;i<prob->l;i++)
        {
                int j = i+rand()%(prob->l-i);
                swap(perm[i],perm[j]);
        }
        for(i=0;i<nr_fold;i++)
        {
                int begin = i*prob->l/nr_fold;
                int end = (i+1)*prob->l/nr_fold;
                int j,k;
                struct svm_problem subprob;

                subprob.l = prob->l-(end-begin);
                subprob.x = Malloc(struct svm_node*,subprob.l);
                subprob.y = Malloc(double,subprob.l);

                k=0;
                for(j=0;j<begin;j++)
                {
                        subprob.x[k] = prob->x[perm[j]];
                        subprob.y[k] = prob->y[perm[j]];
                        ++k;
                }
                for(j=end;j<prob->l;j++)
                {
                        subprob.x[k] = prob->x[perm[j]];
                        subprob.y[k] = prob->y[perm[j]];
                        ++k;
                }
                int p_count=0,n_count=0;
                for(j=0;j<k;j++)
                        if(subprob.y[j]>0)
                                p_count++;
                        else
                                n_count++;

                if(p_count==0 && n_count==0)
                        for(j=begin;j<end;j++)
                                dec_values[perm[j]] = 0;
                else if(p_count > 0 && n_count == 0)
                        for(j=begin;j<end;j++)
                                dec_values[perm[j]] = 1;
                else if(p_count == 0 && n_count > 0)
                        for(j=begin;j<end;j++)
                                dec_values[perm[j]] = -1;
                else
                {
                        svm_parameter subparam = *param;
                        subparam.probability=0;
                        subparam.C=1.0;
                        subparam.nr_weight=2;
                        subparam.weight_label = Malloc(int,2);
                        subparam.weight = Malloc(double,2);
                        subparam.weight_label[0]=+1;
                        subparam.weight_label[1]=-1;
                        subparam.weight[0]=Cp;
                        subparam.weight[1]=Cn;
                        struct svm_model *submodel = svm_train(&subprob,&subparam);
                        for(j=begin;j<end;j++)
                        {
                                svm_predict_values(submodel,prob->x[perm[j]],&(dec_values[perm[j]]));
                                // ensure +1 -1 order; reason not using CV subroutine
                                dec_values[perm[j]] *= submodel->label[0];
                        }
                        svm_destroy_model(submodel);
                        svm_destroy_param(&subparam);
                }
                free(subprob.x);
                free(subprob.y);
        }
        sigmoid_train(prob->l,dec_values,prob->y,probA,probB);
        free(dec_values);
        free(perm);
}

// Return parameter of a Laplace distribution
double svm_svr_probability(
        const svm_problem *prob, const svm_parameter *param)
{
        int i;
        int nr_fold = 5;
        double *ymv = Malloc(double,prob->l);
        double mae = 0;

        svm_parameter newparam = *param;
        newparam.probability = 0;
        svm_cross_validation(prob,&newparam,nr_fold,ymv);
        for(i=0;i<prob->l;i++)
        {
                ymv[i]=prob->y[i]-ymv[i];
                mae += fabs(ymv[i]);
        }
        mae /= prob->l;
        double std=sqrt(2*mae*mae);
        int count=0;
        mae=0;
        for(i=0;i<prob->l;i++)
                if (fabs(ymv[i]) > 5*std)
                        count=count+1;
                else
                        mae+=fabs(ymv[i]);
        mae /= (prob->l-count);
        info("Prob. model for test data: target value = predicted value + z,\nz: Laplace distribution e^(-|z|/sigma)/(2sigma),sigma= %g\n",mae);
        free(ymv);
        return mae;
}


// label: label name, start: begin of each class, count: #data of classes, perm: indices to the original data
// perm, length l, must be allocated before calling this subroutine
void svm_group_classes(const svm_problem *prob, int *nr_class_ret, int **label_ret, int **start_ret, int **count_ret, int *perm)
{
        int l = prob->l;
        int max_nr_class = 16;
        int nr_class = 0;
        int *label = Malloc(int,max_nr_class);
        int *count = Malloc(int,max_nr_class);
        int *data_label = Malloc(int,l);
        int i;

        for(i=0;i<l;i++)
        {
                int this_label = (int)prob->y[i];
                int j;
                for(j=0;j<nr_class;j++)
                {
                        if(this_label == label[j])
                        {
                                ++count[j];
                                break;
                        }
                }
                data_label[i] = j;
                if(j == nr_class)
                {
                        if(nr_class == max_nr_class)
                        {
                                max_nr_class *= 2;
                                label = (int *)realloc(label,max_nr_class*sizeof(int));
                                count = (int *)realloc(count,max_nr_class*sizeof(int));
                        }
                        label[nr_class] = this_label;
                        count[nr_class] = 1;
                        ++nr_class;
                }
        }

        int *start = Malloc(int,nr_class);
        start[0] = 0;
        for(i=1;i<nr_class;i++)
                start[i] = start[i-1]+count[i-1];
        for(i=0;i<l;i++)
        {
                perm[start[data_label[i]]] = i;
                ++start[data_label[i]];
        }
        start[0] = 0;
        for(i=1;i<nr_class;i++)
                start[i] = start[i-1]+count[i-1];

        *nr_class_ret = nr_class;
        *label_ret = label;
        *start_ret = start;
        *count_ret = count;
        free(data_label);
}

//
// Interface functions
//
svm_model *svm_train(const svm_problem *prob, const svm_parameter *param)
{
        svm_model *model = Malloc(svm_model,1);
        model->param = *param;
        model->free_sv = 0;        // XXX

        if(param->svm_type == ONE_CLASS ||
           param->svm_type == EPSILON_SVR ||
           param->svm_type == NU_SVR)
        {
                // regression or one-class-svm
                model->nr_class = 2;
                model->label = NULL;
                model->nSV = NULL;
                model->probA = NULL; model->probB = NULL;
                model->sv_coef = Malloc(double *,1);

                if(param->probability &&
                   (param->svm_type == EPSILON_SVR ||
                    param->svm_type == NU_SVR))
                {
                        model->probA = Malloc(double,1);
                        model->probA[0] = svm_svr_probability(prob,param);
                }

                decision_function f = svm_train_one(prob,param,0,0);
                model->rho = Malloc(double,1);
                model->rho[0] = f.rho;

                int nSV = 0;
                int i;
                for(i=0;i<prob->l;i++)
                        if(fabs(f.alpha[i]) > 0) ++nSV;
                model->l = nSV;
                model->SV = Malloc(svm_node *,nSV);
                model->sv_coef[0] = Malloc(double,nSV);
                int j = 0;
                for(i=0;i<prob->l;i++)
                        if(fabs(f.alpha[i]) > 0)
                        {
                                model->SV[j] = prob->x[i];
                                model->sv_coef[0][j] = f.alpha[i];
                                ++j;
                        }

                free(f.alpha);
        }
        else
        {
                // classification
                int l = prob->l;
                int nr_class;
                int *label = NULL;
                int *start = NULL;
                int *count = NULL;
                int *perm = Malloc(int,l);

                // group training data of the same class
                svm_group_classes(prob,&nr_class,&label,&start,&count,perm);
                svm_node **x = Malloc(svm_node *,l);
                int i;
                for(i=0;i<l;i++)
                        x[i] = prob->x[perm[i]];

                // calculate weighted C

                double *weighted_C = Malloc(double, nr_class);
                for(i=0;i<nr_class;i++)
                        weighted_C[i] = param->C;
                for(i=0;i<param->nr_weight;i++)
                {
                        int j;
                        for(j=0;j<nr_class;j++)
                                if(param->weight_label[i] == label[j])
                                        break;
                        if(j == nr_class)
                                fprintf(stderr,"warning: class label %d specified in weight is not found\n", param->weight_label[i]);
                        else
                                weighted_C[j] *= param->weight[i];
                }

                // train k*(k-1)/2 models

                bool *nonzero = Malloc(bool,l);
                for(i=0;i<l;i++)
                        nonzero[i] = false;
                decision_function *f = Malloc(decision_function,nr_class*(nr_class-1)/2);

                double *probA=NULL,*probB=NULL;
                if (param->probability)
                {
                        probA=Malloc(double,nr_class*(nr_class-1)/2);
                        probB=Malloc(double,nr_class*(nr_class-1)/2);
                }

                int p = 0;
                for(i=0;i<nr_class;i++)
                        for(int j=i+1;j<nr_class;j++)
                        {
                                svm_problem sub_prob;
                                int si = start[i], sj = start[j];
                                int ci = count[i], cj = count[j];
                                sub_prob.l = ci+cj;
                                sub_prob.x = Malloc(svm_node *,sub_prob.l);
                                sub_prob.y = Malloc(double,sub_prob.l);
                                int k;
                                for(k=0;k<ci;k++)
                                {
                                        sub_prob.x[k] = x[si+k];
                                        sub_prob.y[k] = +1;
                                }
                                for(k=0;k<cj;k++)
                                {
                                        sub_prob.x[ci+k] = x[sj+k];
                                        sub_prob.y[ci+k] = -1;
                                }

                                if(param->probability)
                                        svm_binary_svc_probability(&sub_prob,param,weighted_C[i],weighted_C[j],probA[p],probB[p]);

                                f[p] = svm_train_one(&sub_prob,param,weighted_C[i],weighted_C[j]);
                                for(k=0;k<ci;k++)
                                        if(!nonzero[si+k] && fabs(f[p].alpha[k]) > 0)
                                                nonzero[si+k] = true;
                                for(k=0;k<cj;k++)
                                        if(!nonzero[sj+k] && fabs(f[p].alpha[ci+k]) > 0)
                                                nonzero[sj+k] = true;
                                free(sub_prob.x);
                                free(sub_prob.y);
                                ++p;
                        }

                // build output

                model->nr_class = nr_class;

                model->label = Malloc(int,nr_class);
                for(i=0;i<nr_class;i++)
                        model->label[i] = label[i];

                model->rho = Malloc(double,nr_class*(nr_class-1)/2);
                for(i=0;i<nr_class*(nr_class-1)/2;i++)
                        model->rho[i] = f[i].rho;

                if(param->probability)
                {
                        model->probA = Malloc(double,nr_class*(nr_class-1)/2);
                        model->probB = Malloc(double,nr_class*(nr_class-1)/2);
                        for(i=0;i<nr_class*(nr_class-1)/2;i++)
                        {
                                model->probA[i] = probA[i];
                                model->probB[i] = probB[i];
                        }
                }
                else
                {
                        model->probA=NULL;
                        model->probB=NULL;
                }

                int total_sv = 0;
                int *nz_count = Malloc(int,nr_class);
                model->nSV = Malloc(int,nr_class);
                for(i=0;i<nr_class;i++)
                {
                        int nSV = 0;
                        for(int j=0;j<count[i];j++)
                                if(nonzero[start[i]+j])
                                {
                                        ++nSV;
                                        ++total_sv;
                                }
                        model->nSV[i] = nSV;
                        nz_count[i] = nSV;
                }

                info("Total nSV = %d\n",total_sv);

                model->l = total_sv;
                model->SV = Malloc(svm_node *,total_sv);
                p = 0;
                for(i=0;i<l;i++)
                        if(nonzero[i]) model->SV[p++] = x[i];

                int *nz_start = Malloc(int,nr_class);
                nz_start[0] = 0;
                for(i=1;i<nr_class;i++)
                        nz_start[i] = nz_start[i-1]+nz_count[i-1];

                model->sv_coef = Malloc(double *,nr_class-1);
                for(i=0;i<nr_class-1;i++)
                        model->sv_coef[i] = Malloc(double,total_sv);

                p = 0;
                for(i=0;i<nr_class;i++)
                        for(int j=i+1;j<nr_class;j++)
                        {
                                // classifier (i,j): coefficients with
                                // i are in sv_coef[j-1][nz_start[i]...],
                                // j are in sv_coef[i][nz_start[j]...]

                                int si = start[i];
                                int sj = start[j];
                                int ci = count[i];
                                int cj = count[j];

                                int q = nz_start[i];
                                int k;
                                for(k=0;k<ci;k++)
                                        if(nonzero[si+k])
                                                model->sv_coef[j-1][q++] = f[p].alpha[k];
                                q = nz_start[j];
                                for(k=0;k<cj;k++)
                                        if(nonzero[sj+k])
                                                model->sv_coef[i][q++] = f[p].alpha[ci+k];
                                ++p;
                        }

                free(label);
                free(probA);
                free(probB);
                free(count);
                free(perm);
                free(start);
                free(x);
                free(weighted_C);
                free(nonzero);
                for(i=0;i<nr_class*(nr_class-1)/2;i++)
                        free(f[i].alpha);
                free(f);
                free(nz_count);
                free(nz_start);
        }
        return model;
}

// Stratified cross validation
void svm_cross_validation(const svm_problem *prob, const svm_parameter *param, int nr_fold, double *target)
{
        int i;
        int *fold_start = Malloc(int,nr_fold+1);
        int l = prob->l;
        int *perm = Malloc(int,l);
        int nr_class;

        // stratified cv may not give leave-one-out rate
        // Each class to l folds -> some folds may have zero elements
        if((param->svm_type == C_SVC ||
            param->svm_type == NU_SVC) && nr_fold < l)
        {
                int *start = NULL;
                int *label = NULL;
                int *count = NULL;
                svm_group_classes(prob,&nr_class,&label,&start,&count,perm);

                // random shuffle and then data grouped by fold using the array perm
                int *fold_count = Malloc(int,nr_fold);
                int c;
                int *index = Malloc(int,l);
                for(i=0;i<l;i++)
                        index[i]=perm[i];
                for (c=0; c<nr_class; c++)
                        for(i=0;i<count[c];i++)
                        {
                                int j = i+rand()%(count[c]-i);
                                swap(index[start[c]+j],index[start[c]+i]);
                        }
                for(i=0;i<nr_fold;i++)
                {
                        fold_count[i] = 0;
                        for (c=0; c<nr_class;c++)
                                fold_count[i]+=(i+1)*count[c]/nr_fold-i*count[c]/nr_fold;
                }
                fold_start[0]=0;
                for (i=1;i<=nr_fold;i++)
                        fold_start[i] = fold_start[i-1]+fold_count[i-1];
                for (c=0; c<nr_class;c++)
                        for(i=0;i<nr_fold;i++)
                        {
                                int begin = start[c]+i*count[c]/nr_fold;
                                int end = start[c]+(i+1)*count[c]/nr_fold;
                                for(int j=begin;j<end;j++)
                                {
                                        perm[fold_start[i]] = index[j];
                                        fold_start[i]++;
                                }
                        }
                fold_start[0]=0;
                for (i=1;i<=nr_fold;i++)
                        fold_start[i] = fold_start[i-1]+fold_count[i-1];
                free(start);
                free(label);
                free(count);
                free(index);
                free(fold_count);
        }
        else
        {
                for(i=0;i<l;i++) perm[i]=i;
                for(i=0;i<l;i++)
                {
                        int j = i+rand()%(l-i);
                        swap(perm[i],perm[j]);
                }
                for(i=0;i<=nr_fold;i++)
                        fold_start[i]=i*l/nr_fold;
        }

        for(i=0;i<nr_fold;i++)
        {
                int begin = fold_start[i];
                int end = fold_start[i+1];
                int j,k;
                struct svm_problem subprob;

                subprob.l = l-(end-begin);
                subprob.x = Malloc(struct svm_node*,subprob.l);
                subprob.y = Malloc(double,subprob.l);

                k=0;
                for(j=0;j<begin;j++)
                {
                        subprob.x[k] = prob->x[perm[j]];
                        subprob.y[k] = prob->y[perm[j]];
                        ++k;
                }
                for(j=end;j<l;j++)
                {
                        subprob.x[k] = prob->x[perm[j]];
                        subprob.y[k] = prob->y[perm[j]];
                        ++k;
                }
                struct svm_model *submodel = svm_train(&subprob,param);
                if(param->probability &&
                   (param->svm_type == C_SVC || param->svm_type == NU_SVC))
                {
                        double *prob_estimates=Malloc(double,svm_get_nr_class(submodel));
                        for(j=begin;j<end;j++)
                                target[perm[j]] = svm_predict_probability(submodel,prob->x[perm[j]],prob_estimates);
                        free(prob_estimates);
                }
                else
                        for(j=begin;j<end;j++)
                                target[perm[j]] = svm_predict(submodel,prob->x[perm[j]]);
                svm_destroy_model(submodel);
                free(subprob.x);
                free(subprob.y);
        }
        free(fold_start);
        free(perm);
}


int svm_get_svm_type(const svm_model *model)
{
        return model->param.svm_type;
}

int svm_get_nr_class(const svm_model *model)
{
        return model->nr_class;
}

void svm_get_labels(const svm_model *model, int* label)
{
        if (model->label != NULL)
                for(int i=0;i<model->nr_class;i++)
                        label[i] = model->label[i];
}

double svm_get_svr_probability(const svm_model *model)
{
        if ((model->param.svm_type == EPSILON_SVR || model->param.svm_type == NU_SVR) &&
            model->probA!=NULL)
                return model->probA[0];
        else
        {
                info("Model doesn't contain information for SVR probability inference\n");
                return 0;
        }
}

void svm_predict_values(const svm_model *model, const svm_node *x, double* dec_values)
{
        if(model->param.svm_type == ONE_CLASS ||
           model->param.svm_type == EPSILON_SVR ||
           model->param.svm_type == NU_SVR)
        {
                double *sv_coef = model->sv_coef[0];
                double sum = 0;
                for(int i=0;i<model->l;i++)
                        sum += sv_coef[i] * Kernel::k_function(x,model->SV[i],model->param);
                sum -= model->rho[0];
                *dec_values = sum;
        }
        else
        {
                int i;
                int nr_class = model->nr_class;
                int l = model->l;

                double *kvalue = Malloc(double,l);
                for(i=0;i<l;i++)
                        kvalue[i] = Kernel::k_function(x,model->SV[i],model->param);

                int *start = Malloc(int,nr_class);
                start[0] = 0;
                for(i=1;i<nr_class;i++)
                        start[i] = start[i-1]+model->nSV[i-1];

                int p=0;
                for(i=0;i<nr_class;i++)
                        for(int j=i+1;j<nr_class;j++)
                        {
                                double sum = 0;
                                int si = start[i];
                                int sj = start[j];
                                int ci = model->nSV[i];
                                int cj = model->nSV[j];

                                int k;
                                double *coef1 = model->sv_coef[j-1];
                                double *coef2 = model->sv_coef[i];
                                for(k=0;k<ci;k++)
                                        sum += coef1[si+k] * kvalue[si+k];
                                for(k=0;k<cj;k++)
                                        sum += coef2[sj+k] * kvalue[sj+k];
                                sum -= model->rho[p];
                                dec_values[p] = sum;
                                p++;
                        }

                free(kvalue);
                free(start);
        }
}

double svm_predict(const svm_model *model, const svm_node *x)
{
        if(model->param.svm_type == ONE_CLASS ||
           model->param.svm_type == EPSILON_SVR ||
           model->param.svm_type == NU_SVR)
        {
                double res;
                svm_predict_values(model, x, &res);

                if(model->param.svm_type == ONE_CLASS)
                        return (res>0)?1:-1;
                else
                        return res;
        }
        else
        {
                int i;
                int nr_class = model->nr_class;
                double *dec_values = Malloc(double, nr_class*(nr_class-1)/2);
                svm_predict_values(model, x, dec_values);

                int *vote = Malloc(int,nr_class);
                for(i=0;i<nr_class;i++)
                        vote[i] = 0;
                int pos=0;
                for(i=0;i<nr_class;i++)
                        for(int j=i+1;j<nr_class;j++)
                        {
                                if(dec_values[pos++] > 0)
                                        ++vote[i];
                                else
                                        ++vote[j];
                        }

                int vote_max_idx = 0;
                for(i=1;i<nr_class;i++)
                        if(vote[i] > vote[vote_max_idx])
                                vote_max_idx = i;
                free(vote);
                free(dec_values);
                return model->label[vote_max_idx];
        }
}

double svm_predict_probability(
        const svm_model *model, const svm_node *x, double *prob_estimates)
{
        if ((model->param.svm_type == C_SVC || model->param.svm_type == NU_SVC) &&
            model->probA!=NULL && model->probB!=NULL)
        {
                int i;
                int nr_class = model->nr_class;
                double *dec_values = Malloc(double, nr_class*(nr_class-1)/2);
                svm_predict_values(model, x, dec_values);

                double min_prob=1e-7;
                double **pairwise_prob=Malloc(double *,nr_class);
                for(i=0;i<nr_class;i++)
                        pairwise_prob[i]=Malloc(double,nr_class);
                int k=0;
                for(i=0;i<nr_class;i++)
                        for(int j=i+1;j<nr_class;j++)
                        {
                                pairwise_prob[i][j]=min(max(sigmoid_predict(dec_values[k],model->probA[k],model->probB[k]),min_prob),1-min_prob);
                                pairwise_prob[j][i]=1-pairwise_prob[i][j];
                                k++;
                        }
                multiclass_probability(nr_class,pairwise_prob,prob_estimates);

                int prob_max_idx = 0;
                for(i=1;i<nr_class;i++)
                        if(prob_estimates[i] > prob_estimates[prob_max_idx])
                                prob_max_idx = i;
                for(i=0;i<nr_class;i++)
                        free(pairwise_prob[i]);
                free(dec_values);
                free(pairwise_prob);
                return model->label[prob_max_idx];
        }
        else
                return svm_predict(model, x);
}

const char *svm_type_table[] =
{
        "c_svc","nu_svc","one_class","epsilon_svr","nu_svr",NULL
};

const char *kernel_type_table[]=
{
        "linear","polynomial","rbf","sigmoid","precomputed",NULL
};

void ignore_retval_gcc46(int ret) { /* keep gcc46 happy */ }

int svm_save_model(const char *model_file_name, const svm_model *model)
{
        FILE *fp = fopen(model_file_name,"w");
        if(fp==NULL) return -1;

        const svm_parameter& param = model->param;

        fprintf(fp,"svm_type %s\n", svm_type_table[param.svm_type]);
        fprintf(fp,"kernel_type %s\n", kernel_type_table[param.kernel_type]);

        if(param.kernel_type == POLY)
                fprintf(fp,"degree %d\n", param.degree);

        if(param.kernel_type == POLY || param.kernel_type == RBF || param.kernel_type == SIGMOID)
                fprintf(fp,"gamma %g\n", param.gamma);

        if(param.kernel_type == POLY || param.kernel_type == SIGMOID)
                fprintf(fp,"coef0 %g\n", param.coef0);

        int nr_class = model->nr_class;
        int l = model->l;
        fprintf(fp, "nr_class %d\n", nr_class);
        fprintf(fp, "total_sv %d\n",l);

        {
                fprintf(fp, "rho");
                for(int i=0;i<nr_class*(nr_class-1)/2;i++)
                        fprintf(fp," %g",model->rho[i]);
                fprintf(fp, "\n");
        }

        if(model->label)
        {
                fprintf(fp, "label");
                for(int i=0;i<nr_class;i++)
                        fprintf(fp," %d",model->label[i]);
                fprintf(fp, "\n");
        }

        if(model->probA) // regression has probA only
        {
                fprintf(fp, "probA");
                for(int i=0;i<nr_class*(nr_class-1)/2;i++)
                        fprintf(fp," %g",model->probA[i]);
                fprintf(fp, "\n");
        }
        if(model->probB)
        {
                fprintf(fp, "probB");
                for(int i=0;i<nr_class*(nr_class-1)/2;i++)
                        fprintf(fp," %g",model->probB[i]);
                fprintf(fp, "\n");
        }

        if(model->nSV)
        {
                fprintf(fp, "nr_sv");
                for(int i=0;i<nr_class;i++)
                        fprintf(fp," %d",model->nSV[i]);
                fprintf(fp, "\n");
        }

        fprintf(fp, "SV\n");
        const double * const *sv_coef = model->sv_coef;
        const svm_node * const *SV = model->SV;

        for(int i=0;i<l;i++)
        {
                for(int j=0;j<nr_class-1;j++)
                        fprintf(fp, "%.16g ",sv_coef[j][i]);

                const svm_node *p = SV[i];

                if(param.kernel_type == PRECOMPUTED)
                        fprintf(fp,"0:%d ",(int)(p->value));
                else
                        while(p->index != -1)
                        {
                                fprintf(fp,"%d:%.8g ",p->index,p->value);
                                p++;
                        }
                fprintf(fp, "\n");
        }
        if (ferror(fp) != 0 || fclose(fp) != 0) return -1;
        else return 0;
}

svm_model *svm_load_model(const char *model_file_name)
{
        FILE *fp = fopen(model_file_name,"r");
        if(fp==NULL) return NULL;

        // read parameters

        svm_model *model = Malloc(svm_model,1);
        svm_parameter& param = model->param;
        model->rho = NULL;
        model->probA = NULL;
        model->probB = NULL;
        model->label = NULL;
        model->nSV = NULL;

        char cmd[81];
        int ret;
        while(1)
        {
                ret = fscanf(fp,"%80s",cmd);

                if(strcmp(cmd,"svm_type")==0)
                {
                        ret = fscanf(fp,"%80s",cmd);
                        int i;
                        for(i=0;svm_type_table[i];i++)
                        {
                                if(strcmp(svm_type_table[i],cmd)==0)
                                {
                                        param.svm_type=i;
                                        break;
                                }
                        }
                        if(svm_type_table[i] == NULL)
                        {
                                fprintf(stderr,"unknown svm type.\n");
                                free(model->rho);
                                free(model->label);
                                free(model->nSV);
                                free(model);
                                return NULL;
                        }
                }
                else if(strcmp(cmd,"kernel_type")==0)
                {
                        ret = fscanf(fp,"%80s",cmd);
                        int i;
                        for(i=0;kernel_type_table[i];i++)
                        {
                                if(strcmp(kernel_type_table[i],cmd)==0)
                                {
                                        param.kernel_type=i;
                                        break;
                                }
                        }
                        if(kernel_type_table[i] == NULL)
                        {
                                fprintf(stderr,"unknown kernel function.\n");
                                free(model->rho);
                                free(model->label);
                                free(model->nSV);
                                free(model);
                                return NULL;
                        }
                }
                else if(strcmp(cmd,"degree")==0)
                        ret = fscanf(fp,"%d",&param.degree);
                else if(strcmp(cmd,"gamma")==0)
                        ret = fscanf(fp,"%lf",&param.gamma);
                else if(strcmp(cmd,"coef0")==0)
                        ret = fscanf(fp,"%lf",&param.coef0);
                else if(strcmp(cmd,"nr_class")==0)
                        ret = fscanf(fp,"%d",&model->nr_class);
                else if(strcmp(cmd,"total_sv")==0)
                        ret = fscanf(fp,"%d",&model->l);
                else if(strcmp(cmd,"rho")==0)
                {
                        int n = model->nr_class * (model->nr_class-1)/2;
                        model->rho = Malloc(double,n);
                        for(int i=0;i<n;i++)
                                ret = fscanf(fp,"%lf",&model->rho[i]);
                }
                else if(strcmp(cmd,"label")==0)
                {
                        int n = model->nr_class;
                        model->label = Malloc(int,n);
                        for(int i=0;i<n;i++)
                                ret = fscanf(fp,"%d",&model->label[i]);
                }
                else if(strcmp(cmd,"probA")==0)
                {
                        int n = model->nr_class * (model->nr_class-1)/2;
                        model->probA = Malloc(double,n);
                        for(int i=0;i<n;i++)
                                ret = fscanf(fp,"%lf",&model->probA[i]);
                }
                else if(strcmp(cmd,"probB")==0)
                {
                        int n = model->nr_class * (model->nr_class-1)/2;
                        model->probB = Malloc(double,n);
                        for(int i=0;i<n;i++)
                                ret = fscanf(fp,"%lf",&model->probB[i]);
                }
                else if(strcmp(cmd,"nr_sv")==0)
                {
                        int n = model->nr_class;
                        model->nSV = Malloc(int,n);
                        for(int i=0;i<n;i++)
                                ret = fscanf(fp,"%d",&model->nSV[i]);
                }
                else if(strcmp(cmd,"SV")==0)
                {
                        while(1)
                        {
                                int c = getc(fp);
                                if(c==EOF || c=='\n') break;
                        }
                        break;
                }
                else
                {
                        fprintf(stderr,"unknown text in model file: [%s]\n",cmd);
                        free(model->rho);
                        free(model->label);
                        free(model->nSV);
                        free(model);
                        return NULL;
                }
        }

        // read sv_coef and SV

        int elements = 0;
        long pos = ftell(fp);

        while(1)
        {
                int c = fgetc(fp);
                switch(c)
                {
                        case '\n':
                                // count the '-1' element
                        case ':':
                                ++elements;
                                break;
                        case EOF:
                                goto out;
                        default:
                                ;
                }
        }
out:
        fseek(fp,pos,SEEK_SET);

        int m = model->nr_class - 1;
        int l = model->l;
        model->sv_coef = Malloc(double *,m);
        int i;
        for(i=0;i<m;i++)
                model->sv_coef[i] = Malloc(double,l);
        model->SV = Malloc(svm_node*,l);
        svm_node *x_space=NULL;
        if(l>0) x_space = Malloc(svm_node,elements);

        int j=0;
        for(i=0;i<l;i++)
        {
                model->SV[i] = &x_space[j];
                for(int k=0;k<m;k++)
                        ret = fscanf(fp,"%lf",&model->sv_coef[k][i]);
                while(1)
                {
                        int c;
                        do {
                                c = getc(fp);
                                if(c=='\n') goto out2;
                        } while(isspace(c));
                        ungetc(c,fp);
                        ret = fscanf(fp,"%d:%lf",&(x_space[j].index),&(x_space[j].value));
                        ++j;
                }
out2:
                x_space[j++].index = -1;
        }
        if (ferror(fp) != 0 || fclose(fp) != 0) return NULL;

        model->free_sv = 1;        // XXX

        ignore_retval_gcc46(ret);

        return model;
}

void svm_destroy_model(svm_model* model)
{
        if(model->free_sv && model->l > 0)
                free((void *)(model->SV[0]));
        for(int i=0;i<model->nr_class-1;i++)
                free(model->sv_coef[i]);
        free(model->SV);
        free(model->sv_coef);
        free(model->rho);
        free(model->label);
        free(model->probA);
        free(model->probB);
        free(model->nSV);
        free(model);
}

void svm_destroy_param(svm_parameter* param)
{
        free(param->weight_label);
        free(param->weight);
}

const char *svm_check_parameter(const svm_problem *prob, const svm_parameter *param)
{
        // svm_type

        int svm_type = param->svm_type;
        if(svm_type != C_SVC &&
           svm_type != NU_SVC &&
           svm_type != ONE_CLASS &&
           svm_type != EPSILON_SVR &&
           svm_type != NU_SVR)
                return "unknown svm type";

        // kernel_type, degree

        int kernel_type = param->kernel_type;
        if(kernel_type != LINEAR &&
           kernel_type != POLY &&
           kernel_type != RBF &&
           kernel_type != SIGMOID &&
           kernel_type != PRECOMPUTED)
                return "unknown kernel type";

        if(param->degree < 0)
                return "degree of polynomial kernel < 0";

        // cache_size,eps,C,nu,p,shrinking

        if(param->cache_size <= 0)
                return "cache_size <= 0";

        if(param->eps <= 0)
                return "eps <= 0";

        if(svm_type == C_SVC ||
           svm_type == EPSILON_SVR ||
           svm_type == NU_SVR)
                if(param->C <= 0)
                        return "C <= 0";

        if(svm_type == NU_SVC ||
           svm_type == ONE_CLASS ||
           svm_type == NU_SVR)
                if(param->nu <= 0 || param->nu > 1)
                        return "nu <= 0 or nu > 1";

        if(svm_type == EPSILON_SVR)
                if(param->p < 0)
                        return "p < 0";

        if(param->shrinking != 0 &&
           param->shrinking != 1)
                return "shrinking != 0 and shrinking != 1";

        if(param->probability != 0 &&
           param->probability != 1)
                return "probability != 0 and probability != 1";

        if(param->probability == 1 &&
           svm_type == ONE_CLASS)
                return "one-class SVM probability output not supported yet";


        // check whether nu-svc is feasible

        if(svm_type == NU_SVC)
        {
                int l = prob->l;
                int max_nr_class = 16;
                int nr_class = 0;
                int *label = Malloc(int,max_nr_class);
                int *count = Malloc(int,max_nr_class);

                int i;
                for(i=0;i<l;i++)
                {
                        int this_label = (int)prob->y[i];
                        int j;
                        for(j=0;j<nr_class;j++)
                                if(this_label == label[j])
                                {
                                        ++count[j];
                                        break;
                                }
                        if(j == nr_class)
                        {
                                if(nr_class == max_nr_class)
                                {
                                        max_nr_class *= 2;
                                        label = (int *)realloc(label,max_nr_class*sizeof(int));
                                        count = (int *)realloc(count,max_nr_class*sizeof(int));
                                }
                                label[nr_class] = this_label;
                                count[nr_class] = 1;
                                ++nr_class;
                        }
                }

                for(i=0;i<nr_class;i++)
                {
                        int n1 = count[i];
                        for(int j=i+1;j<nr_class;j++)
                        {
                                int n2 = count[j];
                                if(param->nu*(n1+n2)/2 > min(n1,n2))
                                {
                                        free(label);
                                        free(count);
                                        return "specified nu is infeasible";
                                }
                        }
                }
                free(label);
                free(count);
        }

        return NULL;
}

int svm_check_probability_model(const svm_model *model)
{
        return ((model->param.svm_type == C_SVC || model->param.svm_type == NU_SVC) &&
                model->probA!=NULL && model->probB!=NULL) ||
                ((model->param.svm_type == EPSILON_SVR || model->param.svm_type == NU_SVR) &&
                 model->probA!=NULL);
}