svm_c.cpp
上传用户：xgw_05
上传日期：2014-12-08
资源大小：2726k
文件大小：59k
- #include "stdafx.h"
- #include "svm_c.h"
/**
*
* svm_c
*
*/
svm_c::svm_c(){
// initialise arrays
sum =0;
primal=0;
which_alpha = 0;
lambda_eq=0;
sum_alpha = 0;
at_bound=0;
all_alphas=0;
all_ys=0;
working_set=0;
working_set_values=0;
time_init=0;
time_optimize=0;
time_convergence=0;
time_update=0;
time_calc=0;
qp.m = 1;
};
void svm_c::init(kernel_c* new_kernel, parameters_c* new_parameters){
sum =0;
primal=0;
which_alpha = 0;
lambda_eq=0;
at_bound=0;
target_count=0;
kernel = new_kernel;
parameters = new_parameters;
Cpos = parameters->get_Cpos();
Cneg = parameters->get_Cneg();
is_zero = parameters->is_zero;
is_pattern = parameters->is_pattern;
epsilon_pos = parameters->epsilon_pos;
epsilon_neg = parameters->epsilon_neg;
working_set_size = parameters->working_set_size;
sum_alpha = 0;
convergence_epsilon = parameters->convergence_epsilon;
feasible_epsilon = convergence_epsilon;
shrink_const = parameters->shrink_const;
time_init=0;
time_optimize=0;
time_convergence=0;
time_update=0;
time_calc=0;
qp.m = 1;
biased = parameters->biased;
};
void svm_c::init_optimizer(){
if(0 != sum) delete []sum;
if(0 != which_alpha) delete []which_alpha;
if(0 != at_bound) delete []at_bound;
sum = new SVMFLOAT[examples_total];
at_bound = new SVMINT[examples_total];
// init variables
if(working_set_size>examples_total) working_set_size = examples_total;
qp.n = working_set_size;
qp.c = new SVMFLOAT[qp.n];
qp.H = new SVMFLOAT[qp.n*qp.n];
qp.A = new SVMFLOAT[qp.n];
qp.b = new SVMFLOAT[qp.m];
qp.l = new SVMFLOAT[qp.n];
qp.u = new SVMFLOAT[qp.n];
if(! biased){
qp.m = 0;
};
which_alpha = new SVMINT[working_set_size];
primal = new SVMFLOAT[qp.n];
// reserve workspace for calculate_working_set
working_set = new SVMINT[working_set_size];
working_set_values = new SVMFLOAT[working_set_size];
if(parameters->do_scale_y){
epsilon_pos /= examples->get_y_var();
epsilon_neg /= examples->get_y_var();
};
SVMINT i;
// qp.l[i] = 0 done in svm_pattern_c::
for(i=0;i<working_set_size;i++){
qp.l[i] = -is_zero;
};
// Cpos /= (SVMFLOAT)examples_total;
// Cneg /= (SVMFLOAT)examples_total;
if(parameters->quadraticLossPos){
Cpos = infinity;
};
if(parameters->quadraticLossNeg){
Cneg = infinity;
};
// sigfig_max = -log10(is_zero);
lambda_WS = 0;
to_shrink=0;
smo.init(parameters->is_zero,parameters->convergence_epsilon,working_set_size*working_set_size);
};
void svm_c::exit_optimizer(){
delete [](qp.c);
delete [](qp.H);
delete [](qp.A);
delete [](qp.b);
delete [](qp.l);
delete [](qp.u);
delete []primal;
delete []working_set;
delete []working_set_values;
delete []sum;
delete []at_bound;
delete []which_alpha;
primal = 0;
working_set = 0;
working_set_values = 0;
sum=0;
at_bound=0;
which_alpha=0;
};
int svm_c::is_alpha_neg(const SVMINT i){
// variable i is alpha*
// take a look at svm_pattern_c::is_alpha_neg
// and svm_regression_c::is_alpha_neg!
int result=0;
if(is_pattern){
if(all_ys[i] > 0){
result = 1;
}
else{
result = -1;
};
}
else if(all_alphas[i] > 0){
result = 1;
}
else if(all_alphas[i] == 0){
result = 2*(i%2)-1;
}
else{
result = -1;
};
return result;
};
SVMFLOAT svm_c::nabla(const SVMINT i){
if(is_alpha_neg(i) > 0){
return( sum[i] - all_ys[i] + epsilon_neg);
}
else{
return(-sum[i] + all_ys[i] + epsilon_pos);
};
};
SVMFLOAT svm_c::lambda(const SVMINT i){
// size lagrangian multiplier of the active constraint
SVMFLOAT alpha;
SVMFLOAT result = -abs(nabla(i)+is_alpha_neg(i)*lambda_eq);
//= -infinity; // default = not at bound
alpha=all_alphas[i];
if(alpha>is_zero){
// alpha*
if(alpha-Cneg >= - is_zero){
// upper bound active
result = -lambda_eq-nabla(i);
};
}
else if(alpha >= -is_zero){
// lower bound active
if(is_alpha_neg(i) > 0){
result = nabla(i) + lambda_eq;
}
else{
result = nabla(i)-lambda_eq;
};
}
else if(alpha+Cpos <= is_zero){
// upper bound active
result = lambda_eq - nabla(i);
};
return result;
};
int svm_c::feasible(const SVMINT i, SVMFLOAT* the_nabla, SVMFLOAT* the_lambda, int* atbound){
// is direction i feasible to minimize the target function
// (includes which_alpha==0)
int is_feasible=1;
// if(at_bound[i] >= shrink_const){ is_feasible = 0; };
SVMFLOAT alpha;
*the_nabla = nabla(i);
*the_lambda = lambda(i);
alpha=all_alphas[i];
if(alpha-Cneg >= - is_zero){
// alpha* at upper bound
*atbound = 1;
if(*the_lambda >= 0){
at_bound[i]++;
if(at_bound[i] == shrink_const) to_shrink++;
}
else{
at_bound[i] = 0;
};
}
else if((alpha<=is_zero) && (alpha >= -is_zero)){
// lower bound active
*atbound = 1;
if(*the_lambda >= 0){
at_bound[i]++;
if(at_bound[i] == shrink_const) to_shrink++;
}
else{
at_bound[i] = 0;
};
}
else if(alpha+Cpos <= is_zero){
// alpha at upper bound
*atbound = 1;
if(*the_lambda >= 0){
at_bound[i]++;
if(at_bound[i] == shrink_const) to_shrink++;
}
else{
at_bound[i] = 0;
};
}
else{
// not at bound
*atbound = 0;
at_bound[i] = 0;
};
if((*the_lambda >= feasible_epsilon) || (at_bound[i] >= shrink_const)){
is_feasible = 0;
};
return is_feasible;
};
void svm_c::set_svs(example_set_c* training_examples){
// initialised already trained sv (for predicting or testing)
examples = training_examples;
test_set = training_examples;
examples_total = training_examples->size();
all_alphas = examples->get_alphas();
};
void svm_c::reset_shrinked(){
SVMINT old_ex_tot=examples_total;
target_count=0;
examples_total = examples->size();
kernel->set_examples_size(examples_total);
// unshrink, recalculate sum for all variables
SVMINT i,j;
// reset all sums
SVMFLOAT Kij;
for(i=old_ex_tot;i<examples_total;i++){
sum[i] = 0;
at_bound[i] = 0;
};
for(i=0;i<examples_total;i++){
if(abs(all_alphas[i])>is_zero){
for(j=old_ex_tot;j<examples_total;j++){
Kij = kernel->calculate_K(i,j);
sum[j]+=all_alphas[i]*Kij;
};
};
};
sum_alpha=0;
};
SVMFLOAT svm_c::avg_norm2(){
SVMFLOAT avg=0.0;
for(SVMINT i=0;i<examples_total;i++){
avg += kernel->calculate_K(i,i);
};
avg /= (SVMFLOAT)examples_total;
return(avg);
};
svm_result svm_c::train(example_set_c* training_examples){
svm_result the_result;
examples = training_examples;
test_set = training_examples;
if(parameters->verbosity>= 4){
cout<<"training started"<<endl;
};
time_init = get_time();
time_all = time_init;
examples_total = training_examples->size();
all_alphas = examples->get_alphas();
all_ys = examples->get_ys();
SVMINT param_wss = parameters->working_set_size;
if(param_wss > examples_total){
parameters->working_set_size = examples_total;
working_set_size = examples_total;
};
if(parameters->realC <= 0.0){
parameters->realC = 1.0/avg_norm2();
Cpos = parameters->get_Cpos();
Cneg = parameters->get_Cneg();
if(parameters->verbosity>= 2){
cout<<"C set to "<<parameters->realC<<endl;
};
};
if(parameters->balance_cost){
parameters->Lpos = (SVMFLOAT)training_examples->size_pos()/(SVMFLOAT)training_examples->size();
parameters->Lneg = (SVMFLOAT)training_examples->size_neg()/(SVMFLOAT)training_examples->size();
Cpos = parameters->get_Cpos();
Cneg = parameters->get_Cneg();
};
// set up data structure for optimizer
if(parameters->verbosity>= 4){
cout<<"Initing optimizer"<<endl;
};
init_optimizer();
SVMINT iteration = 0;
SVMINT max_iterations = parameters->max_iterations;
if(parameters->verbosity>= 4){
cout<<"Initing working set"<<endl;
};
// WS-optimization
init_working_set();
time_init = get_time() - time_init;
while(iteration < max_iterations){
iteration++;
if(parameters->verbosity>= 5){
cout<<"WS is: ";
// SVMINT uyt;
SVMINT ws_i;
for(ws_i=0;ws_i<working_set_size;ws_i++){
if(is_alpha_neg(working_set[ws_i])>0) cout<<"+"<<working_set[ws_i]<<" ";
if(is_alpha_neg(working_set[ws_i])<0) cout<<"-"<<working_set[ws_i]<<" ";
};
cout<<endl;
cout<<"WS alphas: ";
for(ws_i=0;ws_i<working_set_size;ws_i++){
cout<<all_alphas[working_set[ws_i]];
if((all_alphas[working_set[ws_i]] > Cneg-is_zero) ||
(all_alphas[working_set[ws_i]] < -Cpos+is_zero)){
cout<<"^";
};
cout<<" ";
};
cout<<endl;
cout<<"WS nablas: ";
for(ws_i=0;ws_i<working_set_size;ws_i++){
cout<<nabla(working_set[ws_i])<<" ";
};
cout<<endl;
cout<<"WS lambdas: ";
for(ws_i=0;ws_i<working_set_size;ws_i++){
cout<<lambda(working_set[ws_i])<<" ";
};
cout<<endl;
cout<<"WS at_bounds: ";
for(ws_i=0;ws_i<working_set_size;ws_i++){
cout<<at_bound[working_set[ws_i]]<<" ";
};
cout<<endl;
};
if(parameters->verbosity>= 4){
cout<<"optimizer iteration "<<iteration<<endl;
}
else if(parameters->verbosity>=3){
cout<<".";
cout.flush();
};
optimize();
put_optimizer_values();
SVMINT i;
// SVMFLOAT now_target=0;
// SVMFLOAT now_target_dummy=0;
// for(i=0;i<examples_total;i++){
// now_target_dummy=sum[i]/2-all_ys[i];
// if(is_alpha_neg(i)){
// now_target_dummy+= epsilon_pos;
// }
// else{
// now_target_dummy-= epsilon_neg;
// };
// now_target+=all_alphas[i]*now_target_dummy;
// };
// cout<<"Target function: "<<now_target<<endl;
int conv = convergence();
// old, is_zero is not touched any more
/// @@@
if(conv && (is_zero > parameters->is_zero)){
is_zero = parameters->is_zero;
conv=0;
}
else if(conv) {
if(parameters->verbosity>=3){
// dots
cout<<endl;
};
// check convergence for all parameters
if(0 == is_pattern){
for(i=0;i<examples_total;i++){
if((all_alphas[i]<=is_zero) && (all_alphas[i]>=-is_zero)){
all_alphas[i]=0;
};
};
project_to_constraint();
};
if((examples_total<examples->size()) || (target_count>0)){
// check convergence for all alphas
if(parameters->verbosity>= 2){
cout<<"***** Checking convergence for all variables"<<endl;
};
SVMINT old_target_count=target_count; // t_c set to 0 in reset_shrinked
reset_shrinked();
conv = convergence();
if(0 == conv){
kernel->set_examples_size(examples_total);
target_count = old_target_count;
};
};
if((0 == is_pattern) && (conv)){
// why???
for(i=0;i<examples_total;i++){
if((all_alphas[i]<=is_zero) && (all_alphas[i]>=-is_zero)){
at_bound[i]++;
};
};
conv = convergence();
};
if(conv){
if(parameters->verbosity>= 1){
cout<<"*** Convergence"<<endl;
};
if((parameters->verbosity>=2) ||
(convergence_epsilon > parameters->convergence_epsilon)){
// had to relax the KKT conditions on the way.
// Maybe this isn't necessary any more
SVMFLOAT new_convergence_epsilon = 0;
SVMFLOAT the_lambda;
for(i=0;i<examples_total;i++){
the_lambda = lambda(i);
if(the_lambda < new_convergence_epsilon){
new_convergence_epsilon = the_lambda;
};
};
convergence_epsilon = -new_convergence_epsilon;
};
break;
};
// set variables free again
shrink_const += 10;
for(i=0;i<examples_total;i++){
at_bound[i]=0;
};
};
shrink();
calculate_working_set();
update_working_set();
if(parameters->verbosity >= 4){
SVMINT shrinked=0;
SVMINT upper_bound=0;
SVMINT lower_bound=0;
for(i=0;i<examples_total;i++){
if(at_bound[i] >= shrink_const){
shrinked++;
};
if(abs(all_alphas[i]) < is_zero){
lower_bound++;
}
else if((all_alphas[i]-Cneg>-is_zero) || (all_alphas[i]+Cpos<is_zero)){
upper_bound++;
};
};
cout<<examples_total<<" variables total, ";
cout<<lower_bound<<" variables at lower bound, ";
cout<<upper_bound<<" variables at upper bound, ";
cout<<shrinked<<" variables shrinked"<<endl;
};
};
SVMINT i;
if(iteration >= max_iterations){
cout<<"*** No convergence: Time up."<<endl;
if(examples_total<examples->size()){
// set sums for all variables for statistics
reset_shrinked();
SVMFLOAT new_convergence_epsilon = 0;
SVMFLOAT the_lambda;
for(i=0;i<examples_total;i++){
the_lambda = lambda(i);
if(the_lambda < new_convergence_epsilon){
new_convergence_epsilon = the_lambda;
};
};
convergence_epsilon = -new_convergence_epsilon;
};
};
time_all = get_time() - time_all;
// calculate b
SVMFLOAT new_b=0;
SVMINT new_b_count=0;
for(i=0;i<examples_total;i++){
if((all_alphas[i]-Cneg < -is_zero) && (all_alphas[i]>is_zero)){
new_b += all_ys[i] - sum[i]-epsilon_neg;
new_b_count++;
}
else if((all_alphas[i]+Cpos > is_zero) && (all_alphas[i]<-is_zero)){
new_b += all_ys[i] - sum[i]+epsilon_pos;
new_b_count++;
};
examples->put_alpha(i,all_alphas[i]);
};
if(new_b_count>0){
examples->put_b(new_b/((SVMFLOAT)new_b_count));
}
else{
// unlikely
for(i=0;i<examples_total;i++){
if((all_alphas[i]<is_zero) && (all_alphas[i]>-is_zero)) {
new_b += all_ys[i]- sum[i];
new_b_count++;
};
};
if(new_b_count>0){
examples->put_b(new_b/((SVMFLOAT)new_b_count));
}
else{
// even unlikelier
for(i=0;i<examples_total;i++){
new_b += all_ys[i]- sum[i];
new_b_count++;
};
examples->put_b(new_b/((SVMFLOAT)new_b_count));
};
};
if(parameters->verbosity>= 2){
cout<<"Done training: "<<iteration<<" iterations."<<endl;
if(parameters->verbosity>= 2){
// cout<<"lambda_eq = "<<lambda_eq<<endl;
SVMFLOAT now_target=0;
SVMFLOAT now_target_dummy=0;
for(i=0;i<examples_total;i++){
now_target_dummy=sum[i]/2-all_ys[i];
if(is_alpha_neg(i)){
now_target_dummy+= epsilon_pos;
}
else{
now_target_dummy-= epsilon_neg;
};
now_target+=all_alphas[i]*now_target_dummy;
};
cout<<"Target function: "<<now_target<<endl;
};
};
the_result = print_statistics();
exit_optimizer();
examples->set_initialised_alpha();
parameters->working_set_size = param_wss;
return the_result;
};
void svm_c::shrink(){
// move shrinked examples to back
if(to_shrink>examples_total/10){
SVMINT i;
SVMINT last_pos=examples_total;
if(last_pos > working_set_size){
for(i=0;i<last_pos;i++){
if(at_bound[i] >= shrink_const){
// shrinxk i
sum_alpha += all_alphas[i];
last_pos--;
examples->swap(i,last_pos);
kernel->overwrite(i,last_pos);
sum[i] = sum[last_pos];
at_bound[i] = at_bound[last_pos];
if(last_pos <= working_set_size){
break;
};
};
};
to_shrink=0;
examples_total = last_pos;
kernel->set_examples_size(examples_total);
};
if(parameters->verbosity>=4){
cout<<"shrinked to "<<examples_total<<" variables"<<endl;
};
};
};
int svm_c::convergence(){
long time_start = get_time();
SVMFLOAT the_lambda_eq = 0;
SVMINT total = 0;
SVMFLOAT alpha_sum=0;
SVMFLOAT alpha=0;
SVMINT i;
int result=1;
// actual convergence-test
total = 0; alpha_sum=0;
if(biased){
// calc lambda_eq
for(i=0;i<examples_total;i++){
alpha = all_alphas[i];
alpha_sum += alpha;
if((alpha>is_zero) && (alpha < Cneg)){ //-Cneg < -is_zero)){
// alpha^* = - nabla
the_lambda_eq += -nabla(i); //all_ys[i]-epsilon_neg-sum[i];
total++;
}
else if((alpha<-is_zero) && (alpha > -Cpos)){ //+Cpos > is_zero)){
// alpha = nabla
the_lambda_eq += nabla(i); //all_ys[i]+epsilon_pos-sum[i];
total++;
};
};
if(parameters->verbosity >= 4){
cout<<"lambda_eq = "<<(the_lambda_eq/total)<<endl;
};
if(total>0){
lambda_eq = the_lambda_eq / total;
}
else{
// keep WS lambda_eq
lambda_eq = lambda_WS; //(lambda_eq+4*lambda_WS)/5;
if(parameters->verbosity>= 4){
cout<<"*** no SVs in convergence(), lambda_eq = "<<lambda_eq<<"."<<endl;
};
};
if(target_count>2){
// estimate lambda from WS
if(target_count>20){
// desperate!
lambda_eq = ((40-target_count)*lambda_eq + (target_count-20)*lambda_WS)/20;
if(parameters->verbosity>=5){
cout<<"Re-Re-calculated lambda from WS: "<<lambda_eq<<endl;
};
if(target_count>40){
// really desperate, kick one example out!
i = working_set[target_count%working_set_size];
if(is_alpha_neg(i) > 0){
lambda_eq = -nabla(i);
}
else{
lambda_eq = nabla(i);
};
if(parameters->verbosity>=5){
cout<<"set lambda_eq to nabla("<<i<<"): "<<lambda_eq<<endl;
};
};
}
else{
lambda_eq = lambda_WS;
if(parameters->verbosity>=5){
cout<<"Re-calculated lambda_eq from WS: "<<lambda_eq<<endl;
};
};
};
// check linear constraint
if(abs(alpha_sum+sum_alpha) > convergence_epsilon){
// equality constraint violated
project_to_constraint();
if(parameters->verbosity>= 4){
cout<<"No convergence: equality constraint violated: |"<<(alpha_sum+sum_alpha)<<"| >> 0"<<endl;
};
result = 0;
};
}
else{
// not biased
lambda_eq = 0.0;
};
i=0;
while(i<examples_total){
if(lambda(i)>=-convergence_epsilon){
i++;
}
else{
result = 0;
break;
};
};
time_convergence += get_time() - time_start;
return result;
};
void svm_c::minheap_heapify(SVMINT start, SVMINT size){
// build heap of array working_set[start:start+size-1]
// (i.e. "size" elements starting at "start"th element)
// minheap = 1 <=> maximal element at root
// (i.e. we build the heap of minimal elements)
// v_a[i] = w_s_v[start-1+i], count beginning at 1
SVMFLOAT* value_array = working_set_values+start-1;
SVMINT* pos_array = working_set+start-1;
int running = 1;
SVMINT pos = 1;
SVMINT left, right, largest;
SVMFLOAT dummyf;
SVMINT dummyi;
while(running){
left = 2*pos;
right = left+1;
if((left<=size) &&
(value_array[left] > value_array[pos]))
largest = left;
else{
largest = pos;
};
if((right<=size) &&
(value_array[right] > value_array[largest])){
largest = right;
};
if(largest == pos){
running = 0;
}
else{
//cout<<"switching "<<pos<<" and "<<largest<<endl;
dummyf = value_array[pos];
dummyi = pos_array[pos];
value_array[pos] = value_array[largest];
pos_array[pos] = pos_array[largest];
value_array[largest] = dummyf;
pos_array[largest] = dummyi;
pos = largest;
};
};
};
void svm_c::maxheap_heapify(SVMINT start, SVMINT size){
// build heap of array working_set[start:start+size-1]
// (i.e. "size" elements starting at "start"th element)
// minheap = 1 <=> maximal element at root
// (i.e. we build the heap of minimal elements)
// v_a[i] = w_s_v[start-1+i], count beginning at 1
SVMFLOAT* value_array = working_set_values+start-1;
SVMINT* pos_array = working_set+start-1;
int running = 1;
SVMINT pos = 1;
SVMINT left, right, largest;
SVMFLOAT dummyf;
SVMINT dummyi;
while(running){
left = 2*pos;
right = left+1;
if((left<=size) &&
(value_array[left] < value_array[pos])){
largest = left;
}
else{
largest = pos;
};
if((right<=size) &&
(value_array[right] < value_array[largest])){
largest = right;
};
if(largest == pos){
running = 0;
}
else{
dummyf = value_array[pos];
dummyi = pos_array[pos];
value_array[pos] = value_array[largest];
pos_array[pos] = pos_array[largest];
value_array[largest] = dummyf;
pos_array[largest] = dummyi;
pos = largest;
};
};
};
SVMINT svm_c::maxheap_add(SVMINT size, const SVMINT element, const SVMFLOAT value){
if(size < (working_set_size/2+working_set_size%2)){
// add to max_heap
working_set_values[working_set_size/2+size] = value;
working_set[working_set_size/2+size] = element;
size++;
if(size == working_set_size/2+working_set_size%2){
// build heap
SVMINT j;
for(j=size;j>0;j--){
maxheap_heapify(working_set_size/2+j-1,size+1-j);
};
};
}
else if(value >= working_set_values[working_set_size/2]){
// replace min of max_heap
working_set_values[working_set_size/2] = value;
working_set[working_set_size/2] = element;
maxheap_heapify(working_set_size/2,size);
};
return size;
};
SVMINT svm_c::minheap_add(SVMINT size, const SVMINT element, const SVMFLOAT value){
if(size<working_set_size/2){
// add to min_heap
working_set_values[size] = value;
working_set[size] = element;
size++;
if(size == working_set_size/2){
// build heap
SVMINT j;
for(j=size;j>0;j--){
minheap_heapify(j-1,size+1-j);
};
};
}
else if(value < working_set_values[0]){
// replace max of min_heap
working_set_values[0] = value;
working_set[0] = element;
minheap_heapify(0,size);
};
return size;
};
void svm_c::calculate_working_set(){
/**
*
* Find top and bottom (w.r.t. in_alpha_neg*nabla) feasible
* variables for working_set
*
*/
long time_start = get_time();
// reset WSS
if(working_set_size < parameters->working_set_size){
working_set_size=parameters->working_set_size;
if(working_set_size>examples_total) working_set_size = examples_total;
};
SVMINT heap_min=0;
SVMINT heap_max=0;
SVMINT i=0;
SVMFLOAT sort_value;
working_set_values[0] = infinity;
working_set_values[working_set_size/2] = -infinity;
SVMFLOAT the_lambda;
SVMFLOAT the_nabla;
int is_feasible;
int atbound;
SVMINT j;
while(i<examples_total) {
is_feasible = feasible(i,&the_nabla,&the_lambda,&atbound);
if(0 != is_feasible){
if(is_alpha_neg(i) > 0){
sort_value = -the_nabla; // - : maximum inconsistency approach
}
else{
sort_value = the_nabla;
};
// add to heaps
heap_min = minheap_add(heap_min,i,sort_value);
heap_max = maxheap_add(heap_max,i,sort_value);
};
i++;
};
if(working_set_values[0] >= working_set_values[working_set_size/2]){
if((heap_min>0) &&
(heap_max>0)){
// there could be the same values in the min- and maxheap,
// sort them out (this is very unlikely)
j=0;
i=0;
while(i<heap_min){
// working_set[i] also in max-heap?
j=working_set_size/2;
while((j<working_set_size/2+heap_max) &&
(working_set[j] != working_set[i])){
j++;
};
if(j<working_set_size/2+heap_max){
// working_set[i] equals working_set[j]
if(heap_min<heap_max){
// remove j from WS
working_set[j] = working_set[working_set_size/2-1+heap_max];
heap_max--;
}
else{
working_set[i] = working_set[heap_min-1];
heap_min--;
};
}
else{
i++;
};
};
};
};
if(heap_min+heap_max < working_set_size) {
// condense WS
for(i=0;i<heap_max;i++){
working_set[heap_min+i] = working_set[working_set_size/2+i];
};
working_set_size = heap_min+heap_max;
};
// if((working_set_size<examples_total) && (working_set_size>0)){
if(target_count>0){
// convergence error on last iteration?
// some more tests on WS
// unlikely to happen, so speed isn't so important
// are all variables at the bound?
SVMINT pos_abs;
int bounded_pos=1;
int bounded_neg=1;
SVMINT pos=0;
while((pos<working_set_size) && ((1 == bounded_pos) || (1 == bounded_neg))){
pos_abs = working_set[pos];
if(is_alpha_neg(pos_abs) > 0){
if(all_alphas[pos_abs]-Cneg < -is_zero){
bounded_pos = 0;
};
if(all_alphas[pos_abs] > is_zero){
bounded_neg = 0;
};
}
else{
if(all_alphas[pos_abs]+Cneg > is_zero){
bounded_neg = 0;
};
if(all_alphas[pos_abs] < -is_zero){
bounded_pos = 0;
};
};
pos++;
};
if(0 != bounded_pos){
// all alphas are at upper bound
// need alpha that can be moved upward
// use alpha with smallest lambda
SVMFLOAT max_lambda = infinity;
SVMINT max_pos=examples_total;
for(pos_abs=0;pos_abs<examples_total;pos_abs++){
if(is_alpha_neg(pos_abs) > 0){
if(all_alphas[pos_abs]-Cneg < -is_zero){
if(lambda(pos_abs) < max_lambda){
max_lambda = lambda(pos_abs);
max_pos = pos_abs;
};
};
}
else{
if(all_alphas[pos_abs] < -is_zero){
if(lambda(pos_abs) < max_lambda){
max_lambda = lambda(pos_abs);
max_pos = pos_abs;
};
};
};
};
if(max_pos<examples_total){
if(working_set_size<parameters->working_set_size){
working_set_size++;
};
working_set[working_set_size-1] = max_pos;
};
}
else if(0 != bounded_neg){
// all alphas are at lower bound
// need alpha that can be moved downward
// use alpha with smallest lambda
SVMFLOAT max_lambda = infinity;
SVMINT max_pos=examples_total;
for(pos_abs=0;pos_abs<examples_total;pos_abs++){
if(is_alpha_neg(pos_abs) > 0){
if(all_alphas[pos_abs] > is_zero){
if(lambda(pos_abs) < max_lambda){
max_lambda = lambda(pos_abs);
max_pos = pos_abs;
};
};
}
else{
if(all_alphas[pos_abs]+Cneg > is_zero){
if(lambda(pos_abs) < max_lambda){
max_lambda = lambda(pos_abs);
max_pos = pos_abs;
};
};
};
};
if(max_pos<examples_total){
if(working_set_size<parameters->working_set_size){
working_set_size++;
};
working_set[working_set_size-1] = max_pos;
};
};
};
if((working_set_size<parameters->working_set_size) &&
(working_set_size<examples_total)){
// use full working set
SVMINT pos = (SVMINT)((SVMFLOAT)examples_total*rand()/(RAND_MAX+1.0));
int ok;
while((working_set_size<parameters->working_set_size) &&
(working_set_size<examples_total)){
// add pos into WS if it isn't already
ok = 1;
for(i=0;i<working_set_size;i++){
if(working_set[i] == pos){
ok=0;
i = working_set_size;
};
};
if(1 == ok){
working_set[working_set_size] = pos;
working_set_size++;
};
pos = (pos+1)%examples_total;
};
};
SVMINT ipos;
for(ipos=0;ipos<working_set_size;ipos++){
which_alpha[ipos] = is_alpha_neg(working_set[ipos]);
};
time_calc += get_time() - time_start;
return;
};
void svm_c::project_to_constraint(){
// project alphas to match the constraint
if(biased){
SVMFLOAT alpha_sum = sum_alpha;
SVMINT SVcount=0;
SVMFLOAT alpha;
SVMINT i;
for(i=0;i<examples_total;i++){
alpha = all_alphas[i];
alpha_sum += alpha;
if(((alpha>is_zero) && (alpha-Cneg < -is_zero)) ||
((alpha<-is_zero) && (alpha+Cpos > is_zero))){
SVcount++;
};
};
if(SVcount > 0){
// project
alpha_sum /= (SVMFLOAT)SVcount;
for(i=0;i<examples_total;i++){
alpha = all_alphas[i];
if(((alpha>is_zero) && (alpha-Cneg < -is_zero)) ||
((alpha<-is_zero) && (alpha+Cpos > is_zero))){
all_alphas[i] -= alpha_sum;
};
};
};
};
};
void svm_c::init_working_set(){
// calculate sum
SVMINT i,j;
project_to_constraint();
// check bounds!
if(examples->initialised_alpha()){
if(parameters->verbosity >= 2){
cout<<"Initialising variables, this may take some time."<<endl;
};
for(i=0; i<examples_total;i++){
sum[i] = 0;
at_bound[i] = 0;
for(j=0; j<examples_total;j++){
sum[i] += all_alphas[j]*kernel->calculate_K(i,j);
};
};
}
else{
// skip kernel calculation as all alphas = 0
for(i=0; i<examples_total;i++){
sum[i] = 0;
at_bound[i] = 0;
};
};
if(examples->initialised_alpha()){
calculate_working_set();
}
else{
// first working set is random
j=0;
i=0;
while((i<working_set_size) && (j < examples_total)){
working_set[i] = j;
if(is_alpha_neg(j) > 0){
which_alpha[i] = 1;
}
else{
which_alpha[i] = -1;
};
i++;
j++;
};
};
update_working_set();
};
void svm_c::put_optimizer_values(){
// update nabla, sum, examples.
// sum[i] += (primal_j^*-primal_j-alpha_j^*+alpha_j)K(i,j)
// check for |nabla| < is_zero (nabla <-> nabla*)
// cout<<"put_optimizer_values()"<<endl;
SVMINT i=0;
SVMINT j=0;
SVMINT pos_i;
SVMFLOAT the_new_alpha;
SVMFLOAT* kernel_row;
SVMFLOAT alpha_diff;
long time_start = get_time();
pos_i=working_set_size;
while(pos_i>0){
pos_i--;
if(which_alpha[pos_i]>0){
the_new_alpha = primal[pos_i];
}
else{
the_new_alpha = -primal[pos_i];
};
// next three statements: keep this order!
i = working_set[pos_i];
alpha_diff = the_new_alpha-all_alphas[i];
all_alphas[i] = the_new_alpha;
if(alpha_diff != 0){
// update sum ( => nabla)
kernel_row = kernel->get_row(i);
for(j=0;j<examples_total;j++){
sum[j] += alpha_diff*kernel_row[j];
};
};
};
time_update += get_time() - time_start;
};
void svm_c::update_working_set(){
long time_start = get_time();
// setup subproblem
SVMINT i,j;
SVMINT pos_i, pos_j;
SVMFLOAT* kernel_row;
SVMFLOAT sum_WS;
for(pos_i=0;pos_i<working_set_size;pos_i++){
i = working_set[pos_i];
// put row sort_i in hessian
kernel_row = kernel->get_row(i);
sum_WS=0;
// for(pos_j=0;pos_j<working_set_size;pos_j++){
for(pos_j=0;pos_j<pos_i;pos_j++){
j = working_set[pos_j];
// put all elements K(i,j) in hessian, where j in WS
if(((which_alpha[pos_j] < 0) && (which_alpha[pos_i] < 0)) ||
((which_alpha[pos_j] > 0) && (which_alpha[pos_i] > 0))){
// both i and j positive or negative
(qp.H)[pos_i*working_set_size+pos_j] = kernel_row[j];
(qp.H)[pos_j*working_set_size+pos_i] = kernel_row[j];
}
else{
// one of i and j positive, one negative
(qp.H)[pos_i*working_set_size+pos_j] = -kernel_row[j];
(qp.H)[pos_j*working_set_size+pos_i] = -kernel_row[j];
};
};
for(pos_j=0;pos_j<working_set_size;pos_j++){
j = working_set[pos_j];
sum_WS+=all_alphas[j]*kernel_row[j];
};
// set main diagonal
(qp.H)[pos_i*working_set_size+pos_i] = kernel_row[i];
// linear and box constraints
if(which_alpha[pos_i]<0){
// alpha
(qp.A)[pos_i] = -1;
// lin(alpha) = y_i+eps-sum_{i not in WS} alpha_i K_{ij}
// = y_i+eps-sum_i+sum_{i in WS}
(qp.c)[pos_i] = all_ys[i]+epsilon_pos-sum[i]+sum_WS;
primal[pos_i] = -all_alphas[i];
(qp.u)[pos_i] = Cpos;
}
else{
// alpha^*
(qp.A)[pos_i] = 1;
(qp.c)[pos_i] = -all_ys[i]+epsilon_neg+sum[i]-sum_WS;
primal[pos_i] = all_alphas[i];
(qp.u)[pos_i] = Cneg;
};
};
if(parameters->quadraticLossNeg){
for(pos_i=0;pos_i<working_set_size;pos_i++){
if(which_alpha[pos_i]>0){
(qp.H)[pos_i*(working_set_size+1)] += 1/Cneg;
(qp.u)[pos_i] = infinity;
};
};
};
if(parameters->quadraticLossPos){
for(pos_i=0;pos_i<working_set_size;pos_i++){
if(which_alpha[pos_i]<0){
(qp.H)[pos_i*(working_set_size+1)] += 1/Cpos;
(qp.u)[pos_i] = infinity;
};
};
};
time_update += get_time() - time_start;
};
svm_result svm_c::test(example_set_c* test_examples, int verbose){
svm_result the_result;
test_set = test_examples;
SVMINT i;
SVMFLOAT MAE=0;
SVMFLOAT MSE=0;
SVMFLOAT actloss=0;
SVMFLOAT theloss=0;
SVMFLOAT theloss_pos=0;
SVMFLOAT theloss_neg=0;
SVMINT countpos=0;
SVMINT countneg=0;
// for pattern:
SVMINT correct_pos=0;
SVMINT correct_neg=0;
SVMINT total_pos=0;
SVMINT total_neg=0;
SVMFLOAT prediction;
SVMFLOAT y;
svm_example example;
for(i=0;i<test_set->size();i++){
example = test_set->get_example(i);
prediction = predict(example);
y = examples->unscale_y(test_set->get_y(i));
MAE += abs(y-prediction);
MSE += (y-prediction)*(y-prediction);
actloss=loss(prediction,y);
theloss+=actloss;
if(y < prediction-parameters->epsilon_pos){
theloss_pos += actloss;
countpos++;
}
else if(y > prediction+parameters->epsilon_neg){
theloss_neg += actloss;
countneg++;
};
// if pattern!
if(is_pattern){
if(y>0){
if(prediction>0){
correct_pos++;
};
total_pos++;
}
else{
if(prediction<=0){
correct_neg++;
};
total_neg++;
};
};
};
if(countpos != 0){
theloss_pos /= (SVMFLOAT)countpos;
};
if(countneg != 0){
theloss_neg /= (SVMFLOAT)countneg;
};
the_result.MAE = MAE / (SVMFLOAT)test_set->size();
the_result.MSE = MSE / (SVMFLOAT)test_set->size();
the_result.loss = theloss/test_set->size();
the_result.loss_pos = theloss_pos;
the_result.loss_neg = theloss_neg;
the_result.number_svs = 0;
the_result.number_bsv = 0;
if(is_pattern){
the_result.accuracy = ((SVMFLOAT)(correct_pos+correct_neg))/((SVMFLOAT)(total_pos+total_neg));
the_result.precision = ((SVMFLOAT)correct_pos/((SVMFLOAT)(correct_pos+total_neg-correct_neg)));
the_result.recall = ((SVMFLOAT)correct_pos/(SVMFLOAT)total_pos);
}
else{
the_result.accuracy = -1;
the_result.precision = -1;
the_result.recall = -1;
};
if(verbose){
cout << "Average loss : "<<(theloss/test_set->size())<<endl;
cout << "Avg. loss pos : "<<theloss_pos<<"t ("<<countpos<<" occurences)"<<endl;
cout << "Avg. loss neg : "<<theloss_neg<<"t ("<<countneg<<" occurences)"<<endl;
cout << "Mean absolute error : "<<the_result.MAE<<endl;
cout << "Mean squared error : "<<the_result.MSE<<endl;
if(is_pattern){
// output precision, recall and accuracy
cout<<"Accuracy : "<<the_result.accuracy<<endl;
cout<<"Precision : "<<the_result.precision<<endl;
cout<<"Recall : "<<the_result.recall<<endl;
// nice printout ;-)
int rows = (int)(1+log10((SVMFLOAT)(total_pos+total_neg)));
int now_digits = rows+2;
int i,j;
cout<<endl;
cout<<"Predicted values:"<<endl;
cout<<" |";
for(i=0;i<rows;i++){ cout<<" "; };
cout<<"+ |";
for(j=0;j<rows;j++){ cout<<" "; };
cout<<"-"<<endl;
cout<<"---+";
for(i=0;i<now_digits;i++){ cout<<"-"; };
cout<<"-+-";
for(i=0;i<now_digits;i++){ cout<<"-"; };
cout<<endl;
cout<<" + | ";
now_digits=rows-(int)(1+log10((SVMFLOAT)correct_pos))-1;
for(i=0;i<now_digits;i++){ cout<<" "; };
cout<<correct_pos<<" | ";
now_digits=rows-(int)(1+log10((SVMFLOAT)(total_pos-correct_pos)))-1;
for(i=0;i<now_digits;i++){ cout<<" "; };
cout<<total_pos-correct_pos<<" (true pos)"<<endl;
cout<<" - | ";
now_digits=rows-(int)(1+log10((SVMFLOAT)(total_neg-correct_neg)))-1;
for(i=0;i<now_digits;i++){ cout<<" "; };
cout<<(total_neg-correct_neg)<<" | ";
now_digits=rows-(int)(1+log10((SVMFLOAT)correct_neg))-1;
for(i=0;i<now_digits;i++){ cout<<" "; };
cout<<correct_neg<<" (true neg)"<<endl;
cout<<endl;
};
};
return the_result;
};
void svm_c::optimize(){
// optimizer-specific call
// get time
long time_start = get_time();
qp.n = working_set_size;
SVMINT i;
SVMINT j;
// equality constraint
qp.b[0]=0;
for(i=0;i<working_set_size;i++){
qp.b[0] += all_alphas[working_set[i]];
};
// set initial optimization parameters
SVMFLOAT new_target=0;
SVMFLOAT old_target=0;
SVMFLOAT target_tmp;
for(i=0;i<working_set_size;i++){
target_tmp = primal[i]*qp.H[i*working_set_size+i]/2;
for(j=0;j<i;j++){
target_tmp+=primal[j]*qp.H[j*working_set_size+i];
};
target_tmp+=qp.c[i];
old_target+=target_tmp*primal[i];
};
SVMFLOAT new_constraint_sum=0;
SVMFLOAT my_is_zero = is_zero;
SVMINT sv_count=working_set_size;
qp.n = working_set_size;
// optimize
int KKTerror=1;
int convError=0;
smo.set_max_allowed_error(convergence_epsilon);
// loop while some KKT condition is not valid (alpha=0)
int result;
if(biased){
result = smo.smo_solve(&qp,primal);
lambda_WS = smo.get_lambda_eq();
}
else{
result = smo.smo_solve_single(&qp,primal);
lambda_WS = 0.0;
};
/////////// new
SVMINT it=3;
if(! is_pattern){
// iterate optimization 3 times with changed sign on variables, if KKT conditions are not satisfied
SVMFLOAT lambda_lo;
while(KKTerror && (it>0)){
KKTerror = 0;
it--;
for(i=0;i<working_set_size;i++){
if(primal[i]<is_zero){
lambda_lo = epsilon_neg + epsilon_pos - qp.c[i];
for(j=0;j<working_set_size;j++){
lambda_lo -= primal[j]*qp.H[i*working_set_size+j];
};
if(qp.A[i] > 0){
lambda_lo -= lambda_WS;
}
else{
lambda_lo += lambda_WS;
};
if(lambda_lo<-convergence_epsilon){
// change sign of i
KKTerror=1;
qp.A[i] = -qp.A[i];
which_alpha[i] = -which_alpha[i];
primal[i] = -primal[i];
qp.c[i] = epsilon_neg + epsilon_pos - qp.c[i];
if(qp.A[i]>0){
qp.u[i] = Cneg;
}
else{
qp.u[i] = Cpos;
};
for(j=0;j<working_set_size;j++){
qp.H[i*working_set_size+j] = -qp.H[i*working_set_size+j];
qp.H[j*working_set_size+i] = -qp.H[j*working_set_size+i];
};
if(parameters->quadraticLossNeg){
if(which_alpha[i]>0){
(qp.H)[i*(working_set_size+1)] += 1/Cneg;
(qp.u)[i] = infinity;
}
else{
// previous was neg
(qp.H)[i*(working_set_size+1)] -= 1/Cneg;
};
};
if(parameters->quadraticLossPos){
if(which_alpha[i]<0){
(qp.H)[i*(working_set_size+1)] += 1/Cpos;
(qp.u)[i] = infinity;
}
else{
//previous was pos
(qp.H)[i*(working_set_size+1)] -= 1/Cpos;
};
};
};
};
};
result = smo.smo_solve(&qp,primal);
if(biased){
lambda_WS = smo.get_lambda_eq();
}; // sonst 0 beibehalten
};
};
KKTerror = 1;
//////////////////////
if(parameters->verbosity>=5){
cout<<"smo ended with result "<<result<<endl;
cout<<"lambda_WS = "<<lambda_WS<<endl;
cout<<"smo: Resulting values:"<<endl;
for(i=0;i<working_set_size;i++){
cout<<i<<": "<<primal[i]<<endl;
};
};
while(KKTerror){
// clip
sv_count=working_set_size;
new_constraint_sum=qp.b[0];
if(biased){
// more checks for feasibility
for(i=0;i<working_set_size;i++){
// check if at bound
if(primal[i] <= my_is_zero){
// at lower bound
primal[i] = qp.l[i];
sv_count--;
}
else if(qp.u[i]-primal[i] <= my_is_zero){
// at upper bound
primal[i] = qp.u[i];
sv_count--;
};
new_constraint_sum -= qp.A[i]*primal[i];
};
// enforce equality constraint
if(sv_count>0){
new_constraint_sum /= (SVMFLOAT)sv_count;
if(parameters->verbosity>=5){
cout<<"adjusting "<<sv_count<<" alphas by "<<new_constraint_sum<<endl;
};
for(i=0;i<working_set_size;i++){
if((primal[i] > qp.l[i]) &&
(primal[i] < qp.u[i])){
// real sv
primal[i] += qp.A[i]*new_constraint_sum;
};
};
}
else if(abs(new_constraint_sum)>(SVMFLOAT)working_set_size*is_zero){
// error, can't get feasible point
if(parameters->verbosity>=5){
cout<<"WARNING: No SVs, constraint_sum = "<<new_constraint_sum<<endl;
};
old_target = -infinity;
//is_ok=0;
convError=1;
};
};
// test descend
new_target=0;
for(i=0;i<working_set_size;i++){
// attention: loqo changes one triangle of H!
target_tmp = primal[i]*qp.H[i*working_set_size+i]/2;
for(j=0;j<i;j++){
target_tmp+=primal[j]*qp.H[j*working_set_size+i];
};
target_tmp+=qp.c[i];
new_target+=target_tmp*primal[i];
};
if(new_target < old_target){
KKTerror = 0;
if(parameters->descend < old_target - new_target){
target_count=0;
}
else{
convError=1;
};
if(parameters->verbosity>=5){
cout<<"descend = "<<old_target-new_target<<endl;
};
}
else if(sv_count > 0){
// less SVs
// set my_is_zero to min_i(primal[i]-qp.l[i], qp.u[i]-primal[i])
my_is_zero = 1e20;
for(i=0;i<working_set_size;i++){
if((primal[i] > qp.l[i]) && (primal[i] < qp.u[i])){
if(primal[i] - qp.l[i] < my_is_zero){
my_is_zero = primal[i]-qp.l[i];
};
if(qp.u[i] - primal[i] < my_is_zero){
my_is_zero = qp.u[i] - primal[i];
};
};
};
if(target_count == 0){
my_is_zero *= 2;
};
if(parameters->verbosity>=5){
cout<<"WARNING: no descend ("<<old_target-new_target
<<" <= "<<parameters->descend
// <<", alpha_diff = "<<alpha_diff
<<"), adjusting is_zero to "<<my_is_zero<<endl;
cout<<"new_target = "<<new_target<<endl;
};
}
else{
// nothing we can do
if(parameters->verbosity>=5){
cout<<"WARNING: no descend ("<<old_target-new_target
<<" <= "<<parameters->descend<<"), stopping."<<endl;
};
KKTerror=0;
convError=1;
};
};
if(1 == convError){
target_count++;
// sigfig_max+=0.05;
if(old_target < new_target){
for(i=0;i<working_set_size;i++){
primal[i] = qp.A[i]*all_alphas[working_set[i]];
};
if(parameters->verbosity>=5){
cout<<"WARNING: Convergence error, restoring old primals"<<endl;
};
};
};
if(target_count>50){
// non-recoverable numerical error
convergence_epsilon*=2;
feasible_epsilon = convergence_epsilon;
// sigfig_max=-log10(is_zero);
if(parameters->verbosity>=1)
cout<<"WARNING: reducing KKT precision to "<<convergence_epsilon<<endl;
target_count=0;
};
time_optimize += get_time() - time_start;
};
void svm_c::predict(example_set_c* test_examples){
test_set = test_examples;
SVMINT i;
SVMFLOAT prediction;
svm_example example;
for(i=0;i<test_set->size();i++){
example = test_set->get_example(i);
prediction = predict(example);
test_set->put_y(i,prediction);
};
test_set->set_initialised_y();
test_set->put_b(examples->get_b());
if(parameters->verbosity>=4){
cout<<"Prediction generated"<<endl;
};
};
SVMFLOAT svm_c::predict(svm_example example){
SVMINT i;
svm_example sv;
SVMFLOAT the_sum=examples->get_b();
for(i=0;i<examples_total;i++){
if(all_alphas[i] != 0){
sv = examples->get_example(i);
the_sum += all_alphas[i]*kernel->calculate_K(sv,example);
};
};
the_sum = examples->unscale_y(the_sum);
if(parameters->use_min_prediction){
if(the_sum < parameters->min_prediction){
the_sum = parameters->min_prediction;
};
};
return the_sum;
};
SVMFLOAT svm_c::predict(SVMINT i){
// return (sum[i]+examples->get_b());
// numerically more stable:
return predict(examples->get_example(i));
};
SVMFLOAT svm_c::loss(SVMINT i){
return loss(predict(i),examples->unscale_y(all_ys[i]));
};
SVMFLOAT svm_c::loss(SVMFLOAT prediction, SVMFLOAT value){
SVMFLOAT theloss = prediction - value;
if(is_pattern){
if(((value > 0) && (prediction > 0)) ||
((value <= 0) && (prediction <= 0))){
theloss = 0;
}
};
if(theloss > parameters->epsilon_pos){
if(parameters->quadraticLossPos){
theloss = parameters->Lpos*(theloss-parameters->epsilon_pos)
*(theloss-parameters->epsilon_pos);
}
else{
theloss = parameters->Lpos*(theloss-parameters->epsilon_pos);
};
}
else if(theloss >= -parameters->epsilon_neg){ theloss = 0; }
else{
if(parameters->quadraticLossNeg){
theloss = parameters->Lneg*(-theloss-parameters->epsilon_neg)
*(-theloss-parameters->epsilon_neg);
}
else{
theloss = parameters->Lneg*(-theloss-parameters->epsilon_neg);
};
};
return theloss;
};
void svm_c::print_special_statistics(){
// nothing special here!
};
svm_result svm_c::print_statistics(){
// # SV, # BSV, pos&neg, Loss, VCdim
// Pattern: Acc, Rec, Pred
if(parameters->verbosity>=2){
cout<<"----------------------------------------"<<endl;
};
svm_result the_result;
if(test_set->size() <= 0){
if(parameters->verbosity>= 0){
cout << "No training set given" << endl;
};
the_result.loss = -1;
return the_result; // undefined
};
SVMINT i;
SVMINT svs = 0;
SVMINT bsv = 0;
SVMFLOAT actloss = 0;
SVMFLOAT theloss = 0;
SVMFLOAT theloss_pos=0;
SVMFLOAT theloss_neg=0;
SVMINT countpos=0;
SVMINT countneg=0;
SVMFLOAT min_alpha=infinity;
SVMFLOAT max_alpha=-infinity;
SVMFLOAT norm_w=0;
SVMFLOAT max_norm_x=0;
SVMFLOAT min_norm_x=1e20;
SVMFLOAT norm_x=0;
SVMFLOAT loo_loss_estim=0;
// for pattern:
SVMFLOAT correct_pos=0;
SVMINT correct_neg=0;
SVMINT total_pos=0;
SVMINT total_neg=0;
SVMINT estim_pos=0;
SVMINT estim_neg=0;
SVMFLOAT MSE=0;
SVMFLOAT MAE=0;
SVMFLOAT alpha;
SVMFLOAT prediction;
SVMFLOAT y;
SVMFLOAT xi;
for(i=0;i<examples_total;i++){
// needed before test-loop for performance estimators
norm_w+=all_alphas[i]*sum[i];
alpha=all_alphas[i];
if(alpha!=0){
norm_x = kernel->calculate_K(i,i);
if(norm_x>max_norm_x){
max_norm_x = norm_x;
};
if(norm_x<min_norm_x){
min_norm_x = norm_x;
};
};
};
SVMFLOAT r_delta = max_norm_x;
if(parameters->loo_estim){
r_delta = 0;
SVMFLOAT r_current;
for(SVMINT j=0;j<examples_total;j++){
norm_x = kernel->calculate_K(j,j);
for(i=0;i<examples_total;i++){
r_current = norm_x-kernel->calculate_K(i,j);
if(r_current > r_delta){
r_delta = r_current;
};
};
};
};
for(i=0;i<examples_total;i++){
alpha=all_alphas[i];
if(alpha<min_alpha) min_alpha = alpha;
if(alpha>max_alpha) max_alpha = alpha;
prediction = predict(i);
y = examples->unscale_y(all_ys[i]);
actloss=loss(prediction,y);
theloss+=actloss;
MAE += abs(prediction-y);
MSE += (prediction-y)*(prediction-y);
if(y < prediction-parameters->epsilon_pos){
theloss_pos += actloss;
countpos++;
}
else if(y > prediction+parameters->epsilon_neg){
theloss_neg += actloss;
countneg++;
};
if(parameters->loo_estim){
if(abs(alpha)>is_zero){
if(is_alpha_neg(i)>=0){
loo_loss_estim += loss(prediction-(abs(alpha)*(2*kernel->calculate_K(i,i)+r_delta)+2*epsilon_neg),y);
}
else{
loo_loss_estim += loss(prediction+(abs(alpha)*(2*kernel->calculate_K(i,i)+r_delta)+2*epsilon_pos),y);
};
};
}
else{
// loss doesn't change if non-SV is omitted
loo_loss_estim += actloss;
};
if(abs(alpha)>is_zero){
// a support vector
svs++;
if((alpha-Cneg >= -is_zero) || (alpha+Cpos <= is_zero)){
bsv++;
};
};
if(is_pattern){
if(y>0){
if(prediction>0){
correct_pos++;
};
if(prediction>1){
xi=0;
}
else{
xi=1-prediction;
};
if(2*alpha*r_delta+xi >= 1){
estim_pos++;
};
total_pos++;
}
else{
if(prediction<=0){
correct_neg++;
};
if(prediction<-1){
xi=0;
}
else{
xi=1+prediction;
};
if(2*(-alpha)*r_delta+xi >= 1){
estim_neg++;
};
total_neg++;
};
};
};
if(countpos != 0){
theloss_pos /= (SVMFLOAT)countpos;
};
if(countneg != 0){
theloss_neg /= (SVMFLOAT)countneg;
};
the_result.MAE = MAE / (SVMFLOAT)examples_total;
the_result.MSE = MSE / (SVMFLOAT)examples_total;
the_result.VCdim = 1+norm_w*max_norm_x;
the_result.loss = theloss/((SVMFLOAT)examples_total);
if(parameters->loo_estim){
the_result.pred_loss = loo_loss_estim/((SVMFLOAT)examples_total);
}
else{
the_result.pred_loss = the_result.loss;
};
the_result.loss_pos = theloss_pos;
the_result.loss_neg = theloss_neg;
the_result.number_svs = svs;
the_result.number_bsv = bsv;
if(is_pattern){
the_result.accuracy = ((SVMFLOAT)(correct_pos+correct_neg))/((SVMFLOAT)(total_pos+total_neg));
the_result.precision = ((SVMFLOAT)correct_pos/((SVMFLOAT)(correct_pos+total_neg-correct_neg)));
the_result.recall = ((SVMFLOAT)correct_pos/(SVMFLOAT)total_pos);
if(parameters->loo_estim){
the_result.pred_accuracy = (1-((SVMFLOAT)(estim_pos+estim_neg))/((SVMFLOAT)(total_pos+total_neg)));
the_result.pred_precision = ((SVMFLOAT)(total_pos-estim_pos))/((SVMFLOAT)(total_pos-estim_pos+estim_neg));
the_result.pred_recall = (1-(SVMFLOAT)estim_pos/((SVMFLOAT)total_pos));
}
else{
the_result.pred_accuracy = the_result.accuracy;
the_result.pred_precision = the_result.precision;
the_result.pred_recall = the_result.recall;
};
}
else{
the_result.accuracy = -1;
the_result.precision = -1;
the_result.recall = -1;
the_result.pred_accuracy = -1;
the_result.pred_precision = -1;
the_result.pred_recall = -1;
};
if(convergence_epsilon > parameters->convergence_epsilon){
cout<<"WARNING: The results were obtained using a relaxed epsilon of "<<convergence_epsilon<<" on the KKT conditions!"<<endl;
}
else if(parameters->verbosity>=2){
cout<<"The results are valid with an epsilon of "<<convergence_epsilon<<" on the KKT conditions."<<endl;
};
if(parameters->verbosity >= 2){
cout << "Average loss : "<<the_result.loss<<" (loo-estim: "<< the_result.pred_loss<<")"<<endl;
cout << "Avg. loss pos : "<<theloss_pos<<"t ("<<countpos<<" occurences)"<<endl;
cout << "Avg. loss neg : "<<theloss_neg<<"t ("<<countneg<<" occurences)"<<endl;
cout << "Mean absolute error : "<<the_result.MAE<<endl;
cout << "Mean squared error : "<<the_result.MSE<<endl;
cout << "Support Vectors : "<<svs<<endl;
cout << "Bounded SVs : "<<bsv<<endl;
cout<<"min SV: "<<min_alpha<<endl
<<"max SV: "<<max_alpha<<endl;
cout<<"|w| = "<<sqrt(norm_w)<<endl;
cout<<"max |x| = "<<sqrt(max_norm_x)<<endl;
cout<<"VCdim <= "<<the_result.VCdim<<endl;
print_special_statistics();
if((is_pattern) && (! parameters->is_distribution)){
// output precision, recall and accuracy
if(parameters->loo_estim){
cout<<"performance (+estimators):"<<endl;
cout<<"Accuracy : "<<the_result.accuracy<<" ("<<the_result.pred_accuracy<<")"<<endl;
cout<<"Precision : "<<the_result.precision<<" ("<<the_result.pred_precision<<")"<<endl;
cout<<"Recall : "<<the_result.recall<<" ("<<the_result.pred_recall<<")"<<endl;
}
else{
cout<<"performance :"<<endl;
cout<<"Accuracy : "<<the_result.accuracy<<endl;
cout<<"Precision : "<<the_result.precision<<endl;
cout<<"Recall : "<<the_result.recall<<endl;
};
if(parameters->verbosity>= 2){
// nice printout ;-)
int rows = (int)(1+log10((SVMFLOAT)(total_pos+total_neg)));
int now_digits = rows+2;
int i,j;
cout<<endl;
cout<<"Predicted values:"<<endl;
cout<<" |";
for(i=0;i<rows;i++){ cout<<" "; };
cout<<"+ |";
for(j=0;j<rows;j++){ cout<<" "; };
cout<<"-"<<endl;
cout<<"---+";
for(i=0;i<now_digits;i++){ cout<<"-"; };
cout<<"-+-";
for(i=0;i<now_digits;i++){ cout<<"-"; };
cout<<endl;
cout<<" + | ";
now_digits=rows-(int)(1+log10((SVMFLOAT)correct_pos))-1;
for(i=0;i<now_digits;i++){ cout<<" "; };
cout<<correct_pos<<" | ";
now_digits=rows-(int)(1+log10((SVMFLOAT)(total_pos-correct_pos)))-1;
for(i=0;i<now_digits;i++){ cout<<" "; };
cout<<total_pos-correct_pos<<" (true pos)"<<endl;
cout<<" - | ";
now_digits=rows-(int)(1+log10((SVMFLOAT)(total_neg-correct_neg)))-1;
for(i=0;i<now_digits;i++){ cout<<" "; };
cout<<(total_neg-correct_neg)<<" | ";
now_digits=rows-(int)(1+log10((SVMFLOAT)correct_neg))-1;
for(i=0;i<now_digits;i++){ cout<<" "; };
cout<<correct_neg<<" (true neg)"<<endl;
cout<<endl;
};
};
};
SVMINT dim = examples->get_dim();
if(((parameters->print_w == 1) && (parameters->is_linear != 0)) ||
((dim<100) &&
(parameters->verbosity>= 2) &&
(parameters->is_linear != 0))
){
// print hyperplane
SVMINT j;
svm_example example;
SVMFLOAT* w = new SVMFLOAT[dim];
SVMFLOAT b = examples->get_b();
for(j=0;j<dim;j++) w[j] = 0;
for(i=0;i<examples_total;i++){
example = examples->get_example(i);
alpha = examples->get_alpha(i);
for(j=0;j<example.length;j++){
w[((example.example)[j]).index] += alpha*((example.example)[j]).att;
};
};
if(examples->initialised_scale()){
SVMFLOAT* exp = examples->get_exp();
SVMFLOAT* var = examples->get_var();
for(j=0;j<dim;j++){
if(var[j] != 0){
w[j] /= var[j];
};
if(0 != var[dim]){
w[j] *= var[dim];
};
b -= w[j]*exp[j];
};
b += exp[dim];
};
for(j=0;j<dim;j++){
cout << "w["<<j<<"] = " << w[j] << endl;
};
cout << "b = "<<b<<endl;
if(dim==1){
cout<<"y = "<<w[0]<<"*x+"<<b<<endl;
};
if((dim==2) && (is_pattern)){
cout<<"x1 = "<<-w[0]/w[1]<<"*x0+"<<-b/w[1]<<endl;
};
delete []w;
};
if(parameters->verbosity>= 2){
cout<<"Time for learning:"<<endl
<<"init : "<<(time_init/100)<<"s"<<endl
<<"optimizer : "<<(time_optimize/100)<<"s"<<endl
<<"convergence : "<<(time_convergence/100)<<"s"<<endl
<<"update ws : "<<(time_update/100)<<"s"<<endl
<<"calc ws : "<<(time_calc/100)<<"s"<<endl
<<"============="<<endl
<<"all : "<<(time_all/100)<<"s"<<endl;
}
else if(parameters->verbosity>=2){
cout<<"Time for learning: "<<(time_all/100)<<"s"<<endl;
};
return the_result;
};
/**
*
* Pattern SVM
*
*/
int svm_pattern_c::is_alpha_neg(const SVMINT i){
// variable i is alpha*
int result;
if(all_ys[i] > 0){
result = 1;
}
else{
result = -1;
};
return result;
};
SVMFLOAT svm_pattern_c::nabla(const SVMINT i){
// = is_alpha_neg(i) * sum[i] - 1
if(all_ys[i]>0){
return(sum[i]-1);
}
else{
return(-sum[i]-1);
};
};
SVMFLOAT svm_pattern_c::lambda(const SVMINT i){
// size lagrangian multiplier of the active constraint
SVMFLOAT alpha;
SVMFLOAT result=0;
alpha=all_alphas[i];
if(alpha == Cneg){ //-Cneg >= - is_zero){
// upper bound active
result = -lambda_eq - sum[i] + 1;
}
else if(alpha == 0){ //(alpha <= is_zero) && (alpha >= -is_zero)){
// lower bound active
if(all_ys[i]>0){
result = (sum[i]+lambda_eq) - 1;
}
else{
result = -(sum[i]+lambda_eq) - 1;
};
}
else if(alpha == -Cpos){ //+Cpos <= is_zero){
// upper bound active
result = lambda_eq + sum[i] + 1;
}
else{
if(all_ys[i]>0){
result = -abs(sum[i]+lambda_eq - 1);
}
else{
result = -abs(-sum[i]-lambda_eq - 1);
};
};
return result;
};
int svm_pattern_c::feasible(const SVMINT i, SVMFLOAT* the_nabla, SVMFLOAT* the_lambda, int* atbound){
// is direction i feasible to minimize the target function
// (includes which_alpha==0)
int is_feasible=1;
if(at_bound[i] >= shrink_const){ is_feasible = 0; };
SVMFLOAT alpha;
alpha=all_alphas[i];
if(alpha == Cneg){ //alpha-Cneg >= - is_zero){
// alpha* at upper bound
*atbound = 1;
*the_nabla = sum[i] - 1;
*the_lambda = -lambda_eq - *the_nabla; //sum[i] + 1;
if(*the_lambda >= 0){
at_bound[i]++;
if(at_bound[i] == shrink_const) to_shrink++;
}
else{
at_bound[i] = 0;
};
}
else if(alpha == 0){
// lower bound active
*atbound = -1;
if(all_ys[i]>0){
*the_nabla = sum[i] - 1;
*the_lambda = lambda_eq + *the_nabla; //sum[i]+lambda_eq - 1;
}
else{
*the_nabla = -sum[i] - 1;
*the_lambda = -lambda_eq + *the_nabla; //-sum[i]-lambda_eq - 1;
};
if(*the_lambda >= 0){
at_bound[i]++;
if(at_bound[i] == shrink_const) to_shrink++;
}
else{
at_bound[i] = 0;
};
}
else if(alpha == -Cpos){ //alpha+Cpos <= is_zero){
*atbound = 1;
*the_nabla = -sum[i] - 1;
*the_lambda = lambda_eq - *the_nabla; //sum[i] - 1;
if(*the_lambda >= 0){
at_bound[i]++;
if(at_bound[i] == shrink_const) to_shrink++;
}
else{
at_bound[i] = 0;
};
}
else{
// not at bound
*atbound = 0;
if(all_ys[i]>0){
*the_nabla = sum[i] - 1;
*the_lambda = -abs(*the_nabla+lambda_eq);
}
else{
*the_nabla = -sum[i] - 1;
*the_lambda = -abs(lambda_eq - *the_nabla);
};
at_bound[i] = 0;
};
if(*the_lambda >= feasible_epsilon){
is_feasible = 0;
};
return is_feasible;
};
void svm_pattern_c::init_optimizer(){
// Cs are dived by examples_total in init_optimizer
svm_c::init_optimizer();
SVMINT i;
for(i=0;i<working_set_size;i++){
qp.l[i] = 0;
};
};
/**
*
* Regression SVM
*
*/
int svm_regression_c::is_alpha_neg(const SVMINT i){
// variable i is alpha*
int result;
if(all_alphas[i] > 0){
result = 1;
}
else if(all_alphas[i] == 0){
if(sum[i] - all_ys[i] + lambda_eq>0){
result = -1;
}
else{
result = 1;
};
// result = 2*((i+at_bound[i])%2)-1;
}
else{
result = -1;
};
return result;
};
SVMFLOAT svm_regression_c::nabla(const SVMINT i){
if(all_alphas[i] > 0){
return( sum[i] - all_ys[i] + epsilon_neg);
}
else if(all_alphas[i] == 0){
if(is_alpha_neg(i)>0){
return( sum[i] - all_ys[i] + epsilon_neg);
}
else{
return(-sum[i] + all_ys[i] + epsilon_pos);
};
}
else{
return(-sum[i] + all_ys[i] + epsilon_pos);
};
};
SVMFLOAT svm_regression_c::lambda(const SVMINT i){
// size lagrangian multiplier of the active constraint
SVMFLOAT alpha;
SVMFLOAT result = -abs(nabla(i)+is_alpha_neg(i)*lambda_eq);
//= -infinity; // default = not at bound
alpha=all_alphas[i];
if(alpha>is_zero){
// alpha*
if(alpha-Cneg >= - is_zero){
// upper bound active
result = -lambda_eq-sum[i] + all_ys[i] - epsilon_neg;
};
}
else if(alpha >= -is_zero){
// lower bound active
if(all_alphas[i] > 0){
result = sum[i] - all_ys[i] + epsilon_neg + lambda_eq;
}
else if(all_alphas[i] == 0){
if(is_alpha_neg(i)>0){
result = sum[i] - all_ys[i] + epsilon_neg + lambda_eq;
}
else{
result = -sum[i] + all_ys[i] + epsilon_pos-lambda_eq;
};
}
else{
result = -sum[i] + all_ys[i] + epsilon_pos-lambda_eq;
};
}
else if(alpha+Cpos <= is_zero){
// upper bound active
result = lambda_eq + sum[i] - all_ys[i] - epsilon_pos;
};
return result;
};
int svm_regression_c::feasible(const SVMINT i, SVMFLOAT* the_nabla, SVMFLOAT* the_lambda, int* atbound){
// is direction i feasible to minimize the target function
// (includes which_alpha==0)
int is_feasible = 1; // <=> at_bound < shrink and lambda < -feas_eps
if(at_bound[i] >= shrink_const){ is_feasible = 0; };
SVMFLOAT alpha;
alpha=all_alphas[i];
if(alpha-Cneg >= -is_zero){
// alpha* at upper bound
*atbound = 1;
*the_nabla = sum[i] - all_ys[i] + epsilon_neg;
*the_lambda = -lambda_eq- *the_nabla;
if(*the_lambda >= 0){
at_bound[i]++;
if(at_bound[i] == shrink_const) to_shrink++;
}
else{
at_bound[i] = 0;
};
}
else if((alpha<=is_zero) && (alpha >= -is_zero)){
// lower bound active
*atbound = 1;
if(all_alphas[i] > 0){
*the_nabla = sum[i] - all_ys[i] + epsilon_neg;
*the_lambda = *the_nabla + lambda_eq;
}
else if(all_alphas[i]==0){
if(is_alpha_neg(i)>0){
*the_nabla = sum[i] - all_ys[i] + epsilon_neg;
*the_lambda = *the_nabla + lambda_eq;
}
else{
*the_nabla = -sum[i] + all_ys[i] + epsilon_pos;
*the_lambda = *the_nabla - lambda_eq;
};
}
else{
*the_nabla = -sum[i] + all_ys[i] + epsilon_pos;
*the_lambda = *the_nabla - lambda_eq;
};
if(*the_lambda >= 0){
if(all_alphas[i] != 0){
// check both constraints!
all_alphas[i]=0;
if(is_alpha_neg(i)>0){
*the_lambda = *the_nabla + lambda_eq;
}
else{
*the_lambda = *the_nabla - lambda_eq;
};
if(*the_lambda >= 0){
at_bound[i]++;
};
}
else{
at_bound[i]++;
};
if(at_bound[i] == shrink_const) to_shrink++;
}
else{
at_bound[i] = 0;
};
}
else if(alpha+Cpos <= is_zero){
// alpha at upper bound
*atbound = 1;
*the_nabla = -sum[i] + all_ys[i] + epsilon_pos;
*the_lambda = lambda_eq - *the_nabla;
if(*the_lambda >= 0){
at_bound[i]++;
if(at_bound[i] == shrink_const) to_shrink++;
}
else{
at_bound[i] = 0;
};
}
else{
// not at bound
*atbound = 0;
if(is_alpha_neg(i)>0){
*the_nabla = sum[i] - all_ys[i] + epsilon_neg;
*the_lambda = -abs(*the_nabla + lambda_eq);
}
else{
*the_nabla = -sum[i] + all_ys[i] + epsilon_pos;
*the_lambda = -abs(lambda_eq - *the_nabla);
};
};
if(*the_lambda>=feasible_epsilon){
is_feasible = 0;
};
return is_feasible;
};
English
