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