Algorithm-SVM
view release on metacpan or search on metacpan
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;
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];
while(1)
{
fscanf(fp,"%80s",cmd);
if(strcmp(cmd,"svm_type")==0)
{
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)
{
fscanf(fp,"%80s",cmd);
int i;
for(i=0;kernel_type_table[i];i++)
{
int c;
do {
c = getc(fp);
if(c=='\n') goto out2;
} while(isspace(c));
ungetc(c,fp);
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
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])
{
( run in 0.474 second using v1.01-cache-2.11-cpan-13bb782fe5a )