ext/libsvm/svm.cpp in rb-libsvm-1.4.2 vs ext/libsvm/svm.cpp in rb-libsvm-1.4.3

- old
+ new

@@ -306,11 +306,11 @@ { if(px->index > py->index) ++py; else ++px; - } + } } return sum; } double Kernel::k_function(const svm_node *x, const svm_node *y, @@ -335,11 +335,11 @@ ++y; } else { if(x->index > y->index) - { + { sum += y->value * y->value; ++y; } else { @@ -358,19 +358,19 @@ 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 + return 0; // Unreachable } } // An SMO algorithm in Fan et al., JMLR 6(2005), p. 1889--1918 // Solves: @@ -558,11 +558,11 @@ // optimization step int iter = 0; int max_iter = max(10000000, l>INT_MAX/100 ? INT_MAX : 100*l); int counter = min(l,1000)+1; - + while(iter < max_iter) { // show progress and do shrinking if(--counter == 0) @@ -583,15 +583,15 @@ 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); @@ -606,11 +606,11 @@ 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; @@ -688,11 +688,11 @@ // 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; } @@ -788,19 +788,19 @@ // 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(y[t]==+1) { if(!is_upper_bound(t)) if(-G[t] >= Gmax) { Gmax = -G[t]; @@ -872,11 +872,11 @@ } } } } - if(Gmax+Gmax2 < eps) + if(Gmax+Gmax2 < eps || Gmin_idx == -1) return 1; out_i = Gmax_idx; out_j = Gmin_idx; return 0; @@ -893,11 +893,11 @@ } else if(is_lower_bound(i)) { if(y[i]==+1) return(G[i] > Gmax2); - else + else return(G[i] > Gmax1); } else return(false); } @@ -909,39 +909,39 @@ 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(y[i]==+1) { - if(!is_upper_bound(i)) + if(!is_upper_bound(i)) { if(-G[i] >= Gmax1) Gmax1 = -G[i]; } - if(!is_lower_bound(i)) + if(!is_lower_bound(i)) { if(G[i] >= Gmax2) Gmax2 = G[i]; } } - else + else { - if(!is_upper_bound(i)) + if(!is_upper_bound(i)) { if(-G[i] >= Gmax2) Gmax2 = -G[i]; } - if(!is_lower_bound(i)) + if(!is_lower_bound(i)) { if(G[i] >= Gmax1) Gmax1 = G[i]; } } } - if(unshrink == false && Gmax1 + Gmax2 <= eps*10) + if(unshrink == false && Gmax1 + Gmax2 <= eps*10) { unshrink = true; reconstruct_gradient(); active_size = l; info("*"); @@ -1076,11 +1076,11 @@ for(int j=0;j<active_size;j++) { if(y[j]==+1) { - if (!is_lower_bound(j)) + if (!is_lower_bound(j)) { double grad_diff=Gmaxp+G[j]; if (G[j] >= Gmaxp2) Gmaxp2 = G[j]; if (grad_diff > 0) @@ -1124,11 +1124,11 @@ } } } } - if(max(Gmaxp+Gmaxp2,Gmaxn+Gmaxn2) < eps) + if(max(Gmaxp+Gmaxp2,Gmaxn+Gmaxn2) < eps || Gmin_idx == -1) return 1; if (y[Gmin_idx] == +1) out_i = Gmaxp_idx; else @@ -1142,18 +1142,18 @@ { if(is_upper_bound(i)) { if(y[i]==+1) return(-G[i] > Gmax1); - else + else return(-G[i] > Gmax4); } else if(is_lower_bound(i)) { if(y[i]==+1) return(G[i] > Gmax2); - else + else return(G[i] > Gmax3); } else return(false); } @@ -1178,18 +1178,18 @@ 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]; } } - if(unshrink == false && max(Gmax1+Gmax2,Gmax3+Gmax4) <= eps*10) + if(unshrink == false && max(Gmax1+Gmax2,Gmax3+Gmax4) <= eps*10) { unshrink = true; reconstruct_gradient(); active_size = l; } @@ -1248,36 +1248,36 @@ 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 double[prob.l]; for(int i=0;i<prob.l;i++) QD[i] = (this->*kernel_function)(i,i); } - + Qfloat *get_Q(int i, int len) const { Qfloat *data; int start, j; if((start = cache->get_data(i,&data,len)) < len) @@ -1322,11 +1322,11 @@ cache = new Cache(prob.l,(long int)(param.cache_size*(1<<20))); QD = new double[prob.l]; for(int i=0;i<prob.l;i++) QD[i] = (this->*kernel_function)(i,i); } - + Qfloat *get_Q(int i, int len) const { Qfloat *data; int start, j; if((start = cache->get_data(i,&data,len)) < len) @@ -1358,11 +1358,11 @@ Cache *cache; double *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; @@ -1388,11 +1388,11 @@ { 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 j, real_i = index[i]; if(cache->get_data(real_i,&data,l) < l) @@ -1701,31 +1701,31 @@ return f; } // Platt's binary SVM Probablistic Output: an improvement from Lin et al. static void sigmoid_train( - int l, const double *dec_values, const double *labels, + 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++) @@ -1831,11 +1831,11 @@ 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; @@ -1867,11 +1867,11 @@ 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); @@ -1914,11 +1914,11 @@ 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]]; @@ -1962,23 +1962,23 @@ 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_free_and_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 +// Return parameter of a Laplace distribution static double svm_svr_probability( const svm_problem *prob, const svm_parameter *param) { int i; int nr_fold = 5; @@ -1990,19 +1990,19 @@ 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) + if (fabs(ymv[i]) > 5*std) count=count+1; - else + 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; @@ -2047,12 +2047,12 @@ ++nr_class; } } // - // Labels are ordered by their first occurrence in the training set. - // However, for two-class sets with -1/+1 labels and -1 appears first, + // Labels are ordered by their first occurrence in the training set. + // However, for two-class sets with -1/+1 labels and -1 appears first, // we swap labels to ensure that internally the binary SVM has positive data corresponding to the +1 instances. // if (nr_class == 2 && label[0] == -1 && label[1] == 1) { swap(label[0],label[1]); @@ -2104,11 +2104,11 @@ model->label = NULL; model->nSV = NULL; model->probA = NULL; model->probB = NULL; model->sv_coef = Malloc(double *,1); - if(param->probability && + 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); @@ -2132,11 +2132,11 @@ { model->SV[j] = prob->x[i]; model->sv_coef[0][j] = f.alpha[i]; model->sv_indices[j] = i+1; ++j; - } + } free(f.alpha); } else { @@ -2148,13 +2148,13 @@ 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); - if(nr_class == 1) + if(nr_class == 1) info("WARNING: training data in only one class. See README for details.\n"); - + svm_node **x = Malloc(svm_node *,l); int i; for(i=0;i<l;i++) x[i] = prob->x[perm[i]]; @@ -2162,11 +2162,11 @@ 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) @@ -2174,11 +2174,11 @@ 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); @@ -2227,15 +2227,15 @@ } // 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) @@ -2260,18 +2260,18 @@ 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); model->sv_indices = Malloc(int,total_sv); @@ -2302,11 +2302,11 @@ 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]; @@ -2314,11 +2314,11 @@ 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); @@ -2363,11 +2363,11 @@ 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 (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]); } @@ -2420,11 +2420,11 @@ 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]]; @@ -2435,11 +2435,11 @@ 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 && + 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); @@ -2449,11 +2449,11 @@ for(j=begin;j<end;j++) target[perm[j]] = svm_predict(submodel,prob->x[perm[j]]); svm_free_and_destroy_model(&submodel); free(subprob.x); free(subprob.y); - } + } free(fold_start); free(perm); } @@ -2519,11 +2519,11 @@ } else { 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); @@ -2542,11 +2542,11 @@ 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]; @@ -2580,11 +2580,11 @@ double *dec_values; if(model->param.svm_type == ONE_CLASS || model->param.svm_type == EPSILON_SVR || model->param.svm_type == NU_SVR) dec_values = Malloc(double, 1); - else + else dec_values = Malloc(double, nr_class*(nr_class-1)/2); double pred_result = svm_predict_values(model, x, dec_values); free(dec_values); return pred_result; } @@ -2610,11 +2610,17 @@ { 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); + if (nr_class == 2) + { + prob_estimates[0] = pairwise_prob[0][1]; + prob_estimates[1] = pairwise_prob[1][0]; + } + else + 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; @@ -2622,11 +2628,11 @@ free(pairwise_prob[i]); free(dec_values); free(pairwise_prob); return model->label[prob_max_idx]; } - else + else return svm_predict(model, x); } static const char *svm_type_table[] = { @@ -2641,11 +2647,14 @@ 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; - char *old_locale = strdup(setlocale(LC_ALL, NULL)); + char *old_locale = setlocale(LC_ALL, NULL); + if (old_locale) { + old_locale = strdup(old_locale); + } setlocale(LC_ALL, "C"); const svm_parameter& param = model->param; fprintf(fp,"svm_type %s\n", svm_type_table[param.svm_type]); @@ -2653,27 +2662,27 @@ 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); + fprintf(fp,"gamma %.17g\n", param.gamma); if(param.kernel_type == POLY || param.kernel_type == SIGMOID) - fprintf(fp,"coef0 %g\n", param.coef0); + fprintf(fp,"coef0 %.17g\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," %.17g",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]); @@ -2682,18 +2691,18 @@ 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," %.17g",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," %.17g",model->probB[i]); fprintf(fp, "\n"); } if(model->nSV) { @@ -2708,11 +2717,11 @@ 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]); + fprintf(fp, "%.17g ",sv_coef[j][i]); const svm_node *p = SV[i]; if(param.kernel_type == PRECOMPUTED) fprintf(fp,"0:%d ",(int)(p->value)); @@ -2762,10 +2771,15 @@ // #define FSCANF(_stream, _format, _var) do{ if (fscanf(_stream, _format, _var) != 1) return false; }while(0) bool read_model_header(FILE *fp, svm_model* model) { svm_parameter& param = model->param; + // parameters for training only won't be assigned, but arrays are assigned as NULL for safety + param.nr_weight = 0; + param.weight_label = NULL; + param.weight = NULL; + char cmd[81]; while(1) { FSCANF(fp,"%80s",cmd); @@ -2786,11 +2800,11 @@ fprintf(stderr,"unknown svm type.\n"); return false; } } else if(strcmp(cmd,"kernel_type")==0) - { + { FSCANF(fp,"%80s",cmd); int i; for(i=0;kernel_type_table[i];i++) { if(strcmp(kernel_type_table[i],cmd)==0) @@ -2799,11 +2813,11 @@ break; } } if(kernel_type_table[i] == NULL) { - fprintf(stderr,"unknown kernel function.\n"); + fprintf(stderr,"unknown kernel function.\n"); return false; } } else if(strcmp(cmd,"degree")==0) FSCANF(fp,"%d",&param.degree); @@ -2873,11 +2887,14 @@ svm_model *svm_load_model(const char *model_file_name) { FILE *fp = fopen(model_file_name,"rb"); if(fp==NULL) return NULL; - char *old_locale = strdup(setlocale(LC_ALL, NULL)); + char *old_locale = setlocale(LC_ALL, NULL); + if (old_locale) { + old_locale = strdup(old_locale); + } setlocale(LC_ALL, "C"); // read parameters svm_model *model = Malloc(svm_model,1); @@ -2885,11 +2902,11 @@ model->probA = NULL; model->probB = NULL; model->sv_indices = NULL; model->label = NULL; model->nSV = NULL; - + // read header if (!read_model_header(fp, model)) { fprintf(stderr, "ERROR: fscanf failed to read model\n"); setlocale(LC_ALL, old_locale); @@ -2898,11 +2915,11 @@ free(model->label); free(model->nSV); free(model); return NULL; } - + // read sv_coef and SV int elements = 0; long pos = ftell(fp); @@ -3035,13 +3052,13 @@ 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 && @@ -3090,11 +3107,11 @@ 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; @@ -3123,10 +3140,10 @@ 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++) {