Changeset 7812:adc9f31166b1 in orange


Ignore:
Timestamp:
04/06/11 11:22:42 (3 years ago)
Author:
ales_erjavec <ales.erjavec@…>
Branch:
default
Convert:
dbd3141ee24cfeacb10eecb25c941bcceae168fb
Message:

Updated LibLinear version to 1.7 (has 4 new solver types including a multiclass SVM)

Location:
source/orange
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • source/orange/linear.cpp

    r6963 r7812  
    11/* 
    2 Copyright (c) 2007-2008 The LIBLINEAR Project. 
     2Copyright (c) 2007-2010 The LIBLINEAR Project. 
    33All rights reserved. 
    44 
     
    3737Changes: 
    3838    -#include "linear.h" -> #include "linear.ppp" 
     39    - removed #include "tron.h", instead the contents of tron.h are in linear.h 
    3940    -#ifndef block around swap, min and max definition 
    4041*/ 
     
    6162#endif 
    6263#endif 
     64template <class S, class T> static inline void clone(T*& dst, S* src, int n) 
     65{    
     66    dst = new T[n]; 
     67    memcpy((void *)dst,(void *)src,sizeof(T)*n); 
     68} 
    6369 
    6470#define Malloc(type,n) (type *)malloc((n)*sizeof(type)) 
    6571#define INF HUGE_VAL 
    6672 
    67 #if 0 
    68 void info(const char *fmt,...) 
    69 { 
     73 
     74static void print_string_stdout(const char *s) 
     75{ 
     76    fputs(s,stdout); 
     77    fflush(stdout); 
     78} 
     79 
     80static void (*liblinear_print_string) (const char *) = &print_string_stdout; 
     81 
     82#if 1 
     83static void info(const char *fmt,...) 
     84{ 
     85    char buf[BUFSIZ]; 
    7086    va_list ap; 
    7187    va_start(ap,fmt); 
    72     vprintf(fmt,ap); 
     88    vsprintf(buf,fmt,ap); 
    7389    va_end(ap); 
    74 } 
    75 /*void info_flush() 
    76 { 
    77     fflush(stdout); 
    78 }*/ 
    79 extern void info_flush(); 
     90    (*liblinear_print_string)(buf); 
     91} 
    8092#else 
    81 void info(char *fmt,...) {} 
    82 //extern void info(char *fmt,...); 
    83 void info_flush() {} 
    84 //extern void info_flush(); 
     93static void info(const char *fmt,...) {} 
    8594#endif 
    8695 
    87 class l2_lr_fun : public function1 
     96class l2r_lr_fun : public function 
    8897{ 
    8998public: 
    90     l2_lr_fun(const problem *prob, double Cp, double Cn); 
    91     ~l2_lr_fun(); 
    92  
     99    l2r_lr_fun(const problem *prob, double Cp, double Cn); 
     100    ~l2r_lr_fun(); 
     101     
    93102    double fun(double *w); 
    94103    void grad(double *w, double *g); 
    95104    void Hv(double *s, double *Hs); 
    96  
     105     
    97106    int get_nr_variable(void); 
    98  
     107     
    99108private: 
    100109    void Xv(double *v, double *Xv); 
    101110    void XTv(double *v, double *XTv); 
    102  
     111     
    103112    double *C; 
    104113    double *z; 
     
    107116}; 
    108117 
    109 l2_lr_fun::l2_lr_fun(const problem *prob, double Cp, double Cn) 
     118l2r_lr_fun::l2r_lr_fun(const problem *prob, double Cp, double Cn) 
    110119{ 
    111120    int i; 
    112121    int l=prob->l; 
    113122    int *y=prob->y; 
    114  
     123     
    115124    this->prob = prob; 
    116  
     125     
    117126    z = new double[l]; 
    118127    D = new double[l]; 
    119128    C = new double[l]; 
    120  
     129     
    121130    for (i=0; i<l; i++) 
    122131    { 
     
    128137} 
    129138 
    130 l2_lr_fun::~l2_lr_fun() 
     139l2r_lr_fun::~l2r_lr_fun() 
    131140{ 
    132141    delete[] z; 
     
    136145 
    137146 
    138 double l2_lr_fun::fun(double *w) 
     147double l2r_lr_fun::fun(double *w) 
    139148{ 
    140149    int i; 
     
    142151    int *y=prob->y; 
    143152    int l=prob->l; 
    144     int n=prob->n; 
    145  
     153    int w_size=get_nr_variable(); 
     154     
    146155    Xv(w, z); 
    147156    for(i=0;i<l;i++) 
    148157    { 
    149             double yz = y[i]*z[i]; 
     158        double yz = y[i]*z[i]; 
    150159        if (yz >= 0) 
    151                 f += C[i]*log(1 + exp(-yz)); 
     160            f += C[i]*log(1 + exp(-yz)); 
    152161        else 
    153                 f += C[i]*(-yz+log(1 + exp(yz))); 
     162            f += C[i]*(-yz+log(1 + exp(yz))); 
    154163    } 
    155164    f = 2*f; 
    156     for(i=0;i<n;i++) 
     165    for(i=0;i<w_size;i++) 
    157166        f += w[i]*w[i]; 
    158167    f /= 2.0; 
    159  
     168     
    160169    return(f); 
    161170} 
    162171 
    163 void l2_lr_fun::grad(double *w, double *g) 
     172void l2r_lr_fun::grad(double *w, double *g) 
    164173{ 
    165174    int i; 
    166175    int *y=prob->y; 
    167176    int l=prob->l; 
    168     int n=prob->n; 
    169  
     177    int w_size=get_nr_variable(); 
     178     
    170179    for(i=0;i<l;i++) 
    171180    { 
     
    175184    } 
    176185    XTv(z, g); 
    177  
    178     for(i=0;i<n;i++) 
     186     
     187    for(i=0;i<w_size;i++) 
    179188        g[i] = w[i] + g[i]; 
    180189} 
    181190 
    182 int l2_lr_fun::get_nr_variable(void) 
     191int l2r_lr_fun::get_nr_variable(void) 
    183192{ 
    184193    return prob->n; 
    185194} 
    186195 
    187 void l2_lr_fun::Hv(double *s, double *Hs) 
     196void l2r_lr_fun::Hv(double *s, double *Hs) 
    188197{ 
    189198    int i; 
    190199    int l=prob->l; 
    191     int n=prob->n; 
     200    int w_size=get_nr_variable(); 
    192201    double *wa = new double[l]; 
    193  
     202     
    194203    Xv(s, wa); 
    195204    for(i=0;i<l;i++) 
    196205        wa[i] = C[i]*D[i]*wa[i]; 
    197  
     206     
    198207    XTv(wa, Hs); 
    199     for(i=0;i<n;i++) 
     208    for(i=0;i<w_size;i++) 
    200209        Hs[i] = s[i] + Hs[i]; 
    201210    delete[] wa; 
    202211} 
    203212 
    204 void l2_lr_fun::Xv(double *v, double *Xv) 
     213void l2r_lr_fun::Xv(double *v, double *Xv) 
    205214{ 
    206215    int i; 
    207216    int l=prob->l; 
    208217    feature_node **x=prob->x; 
    209  
     218     
    210219    for(i=0;i<l;i++) 
    211220    { 
     
    220229} 
    221230 
    222 void l2_lr_fun::XTv(double *v, double *XTv) 
     231void l2r_lr_fun::XTv(double *v, double *XTv) 
    223232{ 
    224233    int i; 
    225234    int l=prob->l; 
    226     int n=prob->n; 
     235    int w_size=get_nr_variable(); 
    227236    feature_node **x=prob->x; 
    228  
    229     for(i=0;i<n;i++) 
     237     
     238    for(i=0;i<w_size;i++) 
    230239        XTv[i]=0; 
    231240    for(i=0;i<l;i++) 
     
    240249} 
    241250 
    242 class l2loss_svm_fun : public function1 
     251class l2r_l2_svc_fun : public function 
    243252{ 
    244253public: 
    245     l2loss_svm_fun(const problem *prob, double Cp, double Cn); 
    246     ~l2loss_svm_fun(); 
    247  
     254    l2r_l2_svc_fun(const problem *prob, double Cp, double Cn); 
     255    ~l2r_l2_svc_fun(); 
     256     
    248257    double fun(double *w); 
    249258    void grad(double *w, double *g); 
    250259    void Hv(double *s, double *Hs); 
    251  
     260     
    252261    int get_nr_variable(void); 
    253  
     262     
    254263private: 
    255264    void Xv(double *v, double *Xv); 
    256265    void subXv(double *v, double *Xv); 
    257266    void subXTv(double *v, double *XTv); 
    258  
     267     
    259268    double *C; 
    260269    double *z; 
     
    265274}; 
    266275 
    267 l2loss_svm_fun::l2loss_svm_fun(const problem *prob, double Cp, double Cn) 
     276l2r_l2_svc_fun::l2r_l2_svc_fun(const problem *prob, double Cp, double Cn) 
    268277{ 
    269278    int i; 
    270279    int l=prob->l; 
    271280    int *y=prob->y; 
    272  
     281     
    273282    this->prob = prob; 
    274  
     283     
    275284    z = new double[l]; 
    276285    D = new double[l]; 
    277286    C = new double[l]; 
    278287    I = new int[l]; 
    279  
     288     
    280289    for (i=0; i<l; i++) 
    281290    { 
     
    287296} 
    288297 
    289 l2loss_svm_fun::~l2loss_svm_fun() 
     298l2r_l2_svc_fun::~l2r_l2_svc_fun() 
    290299{ 
    291300    delete[] z; 
     
    295304} 
    296305 
    297 double l2loss_svm_fun::fun(double *w) 
     306double l2r_l2_svc_fun::fun(double *w) 
    298307{ 
    299308    int i; 
     
    301310    int *y=prob->y; 
    302311    int l=prob->l; 
    303     int n=prob->n; 
    304  
     312    int w_size=get_nr_variable(); 
     313     
    305314    Xv(w, z); 
    306315    for(i=0;i<l;i++) 
    307316    { 
    308             z[i] = y[i]*z[i]; 
    309         double d = z[i]-1; 
    310         if (d < 0) 
     317        z[i] = y[i]*z[i]; 
     318        double d = 1-z[i]; 
     319        if (d > 0) 
    311320            f += C[i]*d*d; 
    312321    } 
    313322    f = 2*f; 
    314     for(i=0;i<n;i++) 
     323    for(i=0;i<w_size;i++) 
    315324        f += w[i]*w[i]; 
    316325    f /= 2.0; 
    317  
     326     
    318327    return(f); 
    319328} 
    320329 
    321 void l2loss_svm_fun::grad(double *w, double *g) 
     330void l2r_l2_svc_fun::grad(double *w, double *g) 
    322331{ 
    323332    int i; 
    324333    int *y=prob->y; 
    325334    int l=prob->l; 
    326     int n=prob->n; 
    327  
     335    int w_size=get_nr_variable(); 
     336     
    328337    sizeI = 0; 
    329338    for (i=0;i<l;i++) 
     
    335344        } 
    336345    subXTv(z, g); 
    337  
    338     for(i=0;i<n;i++) 
     346     
     347    for(i=0;i<w_size;i++) 
    339348        g[i] = w[i] + 2*g[i]; 
    340349} 
    341350 
    342 int l2loss_svm_fun::get_nr_variable(void) 
     351int l2r_l2_svc_fun::get_nr_variable(void) 
    343352{ 
    344353    return prob->n; 
    345354} 
    346355 
    347 void l2loss_svm_fun::Hv(double *s, double *Hs) 
     356void l2r_l2_svc_fun::Hv(double *s, double *Hs) 
    348357{ 
    349358    int i; 
    350359    int l=prob->l; 
    351     int n=prob->n; 
     360    int w_size=get_nr_variable(); 
    352361    double *wa = new double[l]; 
    353  
     362     
    354363    subXv(s, wa); 
    355364    for(i=0;i<sizeI;i++) 
    356365        wa[i] = C[I[i]]*wa[i]; 
    357  
     366     
    358367    subXTv(wa, Hs); 
    359     for(i=0;i<n;i++) 
     368    for(i=0;i<w_size;i++) 
    360369        Hs[i] = s[i] + 2*Hs[i]; 
    361370    delete[] wa; 
    362371} 
    363372 
    364 void l2loss_svm_fun::Xv(double *v, double *Xv) 
     373void l2r_l2_svc_fun::Xv(double *v, double *Xv) 
    365374{ 
    366375    int i; 
    367376    int l=prob->l; 
    368377    feature_node **x=prob->x; 
    369  
     378     
    370379    for(i=0;i<l;i++) 
    371380    { 
     
    380389} 
    381390 
    382 void l2loss_svm_fun::subXv(double *v, double *Xv) 
     391void l2r_l2_svc_fun::subXv(double *v, double *Xv) 
    383392{ 
    384393    int i; 
    385394    feature_node **x=prob->x; 
    386  
     395     
    387396    for(i=0;i<sizeI;i++) 
    388397    { 
     
    397406} 
    398407 
    399 void l2loss_svm_fun::subXTv(double *v, double *XTv) 
     408void l2r_l2_svc_fun::subXTv(double *v, double *XTv) 
    400409{ 
    401410    int i; 
    402     int n=prob->n; 
     411    int w_size=get_nr_variable(); 
    403412    feature_node **x=prob->x; 
    404  
    405     for(i=0;i<n;i++) 
     413     
     414    for(i=0;i<w_size;i++) 
    406415        XTv[i]=0; 
    407416    for(i=0;i<sizeI;i++) 
     
    416425} 
    417426 
    418 // A coordinate descent algorithm for 
     427// A coordinate descent algorithm for  
     428// multi-class support vector machines by Crammer and Singer 
     429// 
     430//  min_{\alpha}  0.5 \sum_m ||w_m(\alpha)||^2 + \sum_i \sum_m e^m_i alpha^m_i 
     431//    s.t.     \alpha^m_i <= C^m_i \forall m,i , \sum_m \alpha^m_i=0 \forall i 
     432//  
     433//  where e^m_i = 0 if y_i  = m, 
     434//        e^m_i = 1 if y_i != m, 
     435//  C^m_i = C if m  = y_i,  
     436//  C^m_i = 0 if m != y_i,  
     437//  and w_m(\alpha) = \sum_i \alpha^m_i x_i  
     438// 
     439// Given:  
     440// x, y, C 
     441// eps is the stopping tolerance 
     442// 
     443// solution will be put in w 
     444// 
     445// See Appendix of LIBLINEAR paper, Fan et al. (2008) 
     446 
     447#define GETI(i) (prob->y[i]) 
     448// To support weights for instances, use GETI(i) (i) 
     449 
     450class Solver_MCSVM_CS 
     451{ 
     452public: 
     453    Solver_MCSVM_CS(const problem *prob, int nr_class, double *C, double eps=0.1, int max_iter=100000); 
     454    ~Solver_MCSVM_CS(); 
     455    void Solve(double *w); 
     456private: 
     457    void solve_sub_problem(double A_i, int yi, double C_yi, int active_i, double *alpha_new); 
     458    bool be_shrunk(int i, int m, int yi, double alpha_i, double minG); 
     459    double *B, *C, *G; 
     460    int w_size, l; 
     461    int nr_class; 
     462    int max_iter; 
     463    double eps; 
     464    const problem *prob; 
     465}; 
     466 
     467Solver_MCSVM_CS::Solver_MCSVM_CS(const problem *prob, int nr_class, double *weighted_C, double eps, int max_iter) 
     468{ 
     469    this->w_size = prob->n; 
     470    this->l = prob->l; 
     471    this->nr_class = nr_class; 
     472    this->eps = eps; 
     473    this->max_iter = max_iter; 
     474    this->prob = prob; 
     475    this->B = new double[nr_class]; 
     476    this->G = new double[nr_class]; 
     477    this->C = weighted_C; 
     478} 
     479 
     480Solver_MCSVM_CS::~Solver_MCSVM_CS() 
     481{ 
     482    delete[] B; 
     483    delete[] G; 
     484} 
     485 
     486int compare_double(const void *a, const void *b) 
     487{ 
     488    if(*(double *)a > *(double *)b) 
     489        return -1; 
     490    if(*(double *)a < *(double *)b) 
     491        return 1; 
     492    return 0; 
     493} 
     494 
     495void Solver_MCSVM_CS::solve_sub_problem(double A_i, int yi, double C_yi, int active_i, double *alpha_new) 
     496{ 
     497    int r; 
     498    double *D; 
     499     
     500    clone(D, B, active_i); 
     501    if(yi < active_i) 
     502        D[yi] += A_i*C_yi; 
     503    qsort(D, active_i, sizeof(double), compare_double); 
     504     
     505    double beta = D[0] - A_i*C_yi; 
     506    for(r=1;r<active_i && beta<r*D[r];r++) 
     507        beta += D[r]; 
     508     
     509    beta /= r; 
     510    for(r=0;r<active_i;r++) 
     511    { 
     512        if(r == yi) 
     513            alpha_new[r] = min(C_yi, (beta-B[r])/A_i); 
     514        else 
     515            alpha_new[r] = min((double)0, (beta - B[r])/A_i); 
     516    } 
     517    delete[] D; 
     518} 
     519 
     520bool Solver_MCSVM_CS::be_shrunk(int i, int m, int yi, double alpha_i, double minG) 
     521{ 
     522    double bound = 0; 
     523    if(m == yi) 
     524        bound = C[GETI(i)]; 
     525    if(alpha_i == bound && G[m] < minG) 
     526        return true; 
     527    return false; 
     528} 
     529 
     530void Solver_MCSVM_CS::Solve(double *w) 
     531{ 
     532    int i, m, s; 
     533    int iter = 0; 
     534    double *alpha =  new double[l*nr_class]; 
     535    double *alpha_new = new double[nr_class]; 
     536    int *index = new int[l]; 
     537    double *QD = new double[l]; 
     538    int *d_ind = new int[nr_class]; 
     539    double *d_val = new double[nr_class]; 
     540    int *alpha_index = new int[nr_class*l]; 
     541    int *y_index = new int[l]; 
     542    int active_size = l; 
     543    int *active_size_i = new int[l]; 
     544    double eps_shrink = max(10.0*eps, 1.0); // stopping tolerance for shrinking 
     545    bool start_from_all = true; 
     546    // initial 
     547    for(i=0;i<l*nr_class;i++) 
     548        alpha[i] = 0; 
     549    for(i=0;i<w_size*nr_class;i++) 
     550        w[i] = 0;  
     551    for(i=0;i<l;i++) 
     552    { 
     553        for(m=0;m<nr_class;m++) 
     554            alpha_index[i*nr_class+m] = m; 
     555        feature_node *xi = prob->x[i]; 
     556        QD[i] = 0; 
     557        while(xi->index != -1) 
     558        { 
     559            QD[i] += (xi->value)*(xi->value); 
     560            xi++; 
     561        } 
     562        active_size_i[i] = nr_class; 
     563        y_index[i] = prob->y[i]; 
     564        index[i] = i; 
     565    } 
     566     
     567    while(iter < max_iter)  
     568    { 
     569        double stopping = -INF; 
     570        for(i=0;i<active_size;i++) 
     571        { 
     572            int j = i+rand()%(active_size-i); 
     573            swap(index[i], index[j]); 
     574        } 
     575        for(s=0;s<active_size;s++) 
     576        { 
     577            i = index[s]; 
     578            double Ai = QD[i]; 
     579            double *alpha_i = &alpha[i*nr_class]; 
     580            int *alpha_index_i = &alpha_index[i*nr_class]; 
     581             
     582            if(Ai > 0) 
     583            { 
     584                for(m=0;m<active_size_i[i];m++) 
     585                    G[m] = 1; 
     586                if(y_index[i] < active_size_i[i]) 
     587                    G[y_index[i]] = 0; 
     588                 
     589                feature_node *xi = prob->x[i]; 
     590                while(xi->index!= -1) 
     591                { 
     592                    double *w_i = &w[(xi->index-1)*nr_class]; 
     593                    for(m=0;m<active_size_i[i];m++) 
     594                        G[m] += w_i[alpha_index_i[m]]*(xi->value); 
     595                    xi++; 
     596                } 
     597                 
     598                double minG = INF; 
     599                double maxG = -INF; 
     600                for(m=0;m<active_size_i[i];m++) 
     601                { 
     602                    if(alpha_i[alpha_index_i[m]] < 0 && G[m] < minG) 
     603                        minG = G[m]; 
     604                    if(G[m] > maxG) 
     605                        maxG = G[m]; 
     606                } 
     607                if(y_index[i] < active_size_i[i]) 
     608                    if(alpha_i[prob->y[i]] < C[GETI(i)] && G[y_index[i]] < minG) 
     609                        minG = G[y_index[i]]; 
     610                 
     611                for(m=0;m<active_size_i[i];m++) 
     612                { 
     613                    if(be_shrunk(i, m, y_index[i], alpha_i[alpha_index_i[m]], minG)) 
     614                    { 
     615                        active_size_i[i]--; 
     616                        while(active_size_i[i]>m) 
     617                        { 
     618                            if(!be_shrunk(i, active_size_i[i], y_index[i],  
     619                                          alpha_i[alpha_index_i[active_size_i[i]]], minG)) 
     620                            { 
     621                                swap(alpha_index_i[m], alpha_index_i[active_size_i[i]]); 
     622                                swap(G[m], G[active_size_i[i]]); 
     623                                if(y_index[i] == active_size_i[i]) 
     624                                    y_index[i] = m; 
     625                                else if(y_index[i] == m)  
     626                                    y_index[i] = active_size_i[i]; 
     627                                break; 
     628                            } 
     629                            active_size_i[i]--; 
     630                        } 
     631                    } 
     632                } 
     633                 
     634                if(active_size_i[i] <= 1) 
     635                { 
     636                    active_size--; 
     637                    swap(index[s], index[active_size]); 
     638                    s--;     
     639                    continue; 
     640                } 
     641                 
     642                if(maxG-minG <= 1e-12) 
     643                    continue; 
     644                else 
     645                    stopping = max(maxG - minG, stopping); 
     646                 
     647                for(m=0;m<active_size_i[i];m++) 
     648                    B[m] = G[m] - Ai*alpha_i[alpha_index_i[m]] ; 
     649                 
     650                solve_sub_problem(Ai, y_index[i], C[GETI(i)], active_size_i[i], alpha_new); 
     651                int nz_d = 0; 
     652                for(m=0;m<active_size_i[i];m++) 
     653                { 
     654                    double d = alpha_new[m] - alpha_i[alpha_index_i[m]]; 
     655                    alpha_i[alpha_index_i[m]] = alpha_new[m]; 
     656                    if(fabs(d) >= 1e-12) 
     657                    { 
     658                        d_ind[nz_d] = alpha_index_i[m]; 
     659                        d_val[nz_d] = d; 
     660                        nz_d++; 
     661                    } 
     662                } 
     663                 
     664                xi = prob->x[i]; 
     665                while(xi->index != -1) 
     666                { 
     667                    double *w_i = &w[(xi->index-1)*nr_class]; 
     668                    for(m=0;m<nz_d;m++) 
     669                        w_i[d_ind[m]] += d_val[m]*xi->value; 
     670                    xi++; 
     671                } 
     672            } 
     673        } 
     674         
     675        iter++; 
     676        if(iter % 10 == 0) 
     677        { 
     678            info("."); 
     679        } 
     680         
     681        if(stopping < eps_shrink) 
     682        { 
     683            if(stopping < eps && start_from_all == true) 
     684                break; 
     685            else 
     686            { 
     687                active_size = l; 
     688                for(i=0;i<l;i++) 
     689                    active_size_i[i] = nr_class; 
     690                info("*"); 
     691                eps_shrink = max(eps_shrink/2, eps); 
     692                start_from_all = true; 
     693            } 
     694        } 
     695        else 
     696            start_from_all = false; 
     697    } 
     698     
     699    info("\noptimization finished, #iter = %d\n",iter); 
     700    if (iter >= max_iter) 
     701        info("\nWARNING: reaching max number of iterations\n"); 
     702     
     703    // calculate objective value 
     704    double v = 0; 
     705    int nSV = 0; 
     706    for(i=0;i<w_size*nr_class;i++) 
     707        v += w[i]*w[i]; 
     708    v = 0.5*v; 
     709    for(i=0;i<l*nr_class;i++) 
     710    { 
     711        v += alpha[i]; 
     712        if(fabs(alpha[i]) > 0) 
     713            nSV++; 
     714    } 
     715    for(i=0;i<l;i++) 
     716        v -= alpha[i*nr_class+prob->y[i]]; 
     717    info("Objective value = %lf\n",v); 
     718    info("nSV = %d\n",nSV); 
     719     
     720    delete [] alpha; 
     721    delete [] alpha_new; 
     722    delete [] index; 
     723    delete [] QD; 
     724    delete [] d_ind; 
     725    delete [] d_val; 
     726    delete [] alpha_index; 
     727    delete [] y_index; 
     728    delete [] active_size_i; 
     729} 
     730 
     731// A coordinate descent algorithm for  
    419732// L1-loss and L2-loss SVM dual problems 
    420733// 
    421734//  min_\alpha  0.5(\alpha^T (Q + D)\alpha) - e^T \alpha, 
    422735//    s.t.      0 <= alpha_i <= upper_bound_i, 
    423 // 
     736//  
    424737//  where Qij = yi yj xi^T xj and 
    425 //  D is a diagonal matrix 
     738//  D is a diagonal matrix  
    426739// 
    427740// In L1-SVM case: 
     
    429742//      upper_bound_i = Cn if y_i = -1 
    430743//      D_ii = 0 
    431 // In L2-Svm case: 
     744// In L2-SVM case: 
    432745//      upper_bound_i = INF 
    433746//      D_ii = 1/(2*Cp) if y_i = 1 
    434747//      D_ii = 1/(2*Cn) if y_i = -1 
    435748// 
    436 // Given: 
     749// Given:  
    437750// x, y, Cp, Cn 
    438751// eps is the stopping tolerance 
    439752// 
    440753// solution will be put in w 
    441  
    442 static void solve_linear_c_svc( 
    443     const problem *prob, double *w, double eps, 
    444     double Cp, double Cn, int solver_type) 
     754//  
     755// See Algorithm 3 of Hsieh et al., ICML 2008 
     756 
     757#undef GETI 
     758#define GETI(i) (y[i]+1) 
     759// To support weights for instances, use GETI(i) (i) 
     760 
     761static void solve_l2r_l1l2_svc( 
     762                               const problem *prob, double *w, double eps,  
     763                               double Cp, double Cn, int solver_type) 
    445764{ 
    446765    int l = prob->l; 
    447     int n = prob->n; 
     766    int w_size = prob->n; 
    448767    int i, s, iter = 0; 
    449768    double C, d, G; 
    450769    double *QD = new double[l]; 
    451     int max_iter = 20000; 
     770    int max_iter = 1000; 
    452771    int *index = new int[l]; 
    453772    double *alpha = new double[l]; 
    454773    schar *y = new schar[l]; 
    455774    int active_size = l; 
    456  
     775     
    457776    // PG: projected gradient, for shrinking and stopping 
    458777    double PG; 
     
    460779    double PGmin_old = -INF; 
    461780    double PGmax_new, PGmin_new; 
    462  
    463     // default solver_type: L2LOSS_SVM_DUAL 
    464     double diag_p = 0.5/Cp, diag_n = 0.5/Cn; 
    465     double upper_bound_p = INF, upper_bound_n = INF; 
    466     if(solver_type == L1LOSS_SVM_DUAL) 
    467     { 
    468         diag_p = 0; diag_n = 0; 
    469         upper_bound_p = Cp; upper_bound_n = Cn; 
    470     } 
    471  
    472     for(i=0; i<n; i++) 
     781     
     782    // default solver_type: L2R_L2LOSS_SVC_DUAL 
     783    double diag[3] = {0.5/Cn, 0, 0.5/Cp}; 
     784    double upper_bound[3] = {INF, 0, INF}; 
     785    if(solver_type == L2R_L1LOSS_SVC_DUAL) 
     786    { 
     787        diag[0] = 0; 
     788        diag[2] = 0; 
     789        upper_bound[0] = Cn; 
     790        upper_bound[2] = Cp; 
     791    } 
     792     
     793    for(i=0; i<w_size; i++) 
    473794        w[i] = 0; 
    474795    for(i=0; i<l; i++) 
     
    477798        if(prob->y[i] > 0) 
    478799        { 
    479             y[i] = +1; 
    480             QD[i] = diag_p; 
     800            y[i] = +1;  
    481801        } 
    482802        else 
    483803        { 
    484804            y[i] = -1; 
    485             QD[i] = diag_n; 
    486         } 
    487  
     805        } 
     806        QD[i] = diag[GETI(i)]; 
     807         
    488808        feature_node *xi = prob->x[i]; 
    489809        while (xi->index != -1) 
     
    494814        index[i] = i; 
    495815    } 
    496  
     816     
    497817    while (iter < max_iter) 
    498818    { 
    499819        PGmax_new = -INF; 
    500820        PGmin_new = INF; 
    501  
     821         
    502822        for (i=0; i<active_size; i++) 
    503823        { 
     
    505825            swap(index[i], index[j]); 
    506826        } 
    507  
    508         for (s=0; s < active_size; s++) 
     827         
     828        for (s=0; s<active_size; s++) 
    509829        { 
    510830            i = index[s]; 
    511831            G = 0; 
    512832            schar yi = y[i]; 
    513  
     833             
    514834            feature_node *xi = prob->x[i]; 
    515835            while(xi->index!= -1) 
     
    519839            } 
    520840            G = G*yi-1; 
    521  
    522             if(yi == 1) 
    523             { 
    524                 C = upper_bound_p; 
    525                 G += alpha[i]*diag_p; 
    526             } 
    527             else 
    528             { 
    529                 C = upper_bound_n; 
    530                 G += alpha[i]*diag_n; 
    531             } 
    532  
     841             
     842            C = upper_bound[GETI(i)]; 
     843            G += alpha[i]*diag[GETI(i)]; 
     844             
    533845            PG = 0; 
    534             if (alpha[i] ==0) 
     846            if (alpha[i] == 0) 
    535847            { 
    536848                if (G > PGmax_old) 
     
    558870            else 
    559871                PG = G; 
    560  
     872             
    561873            PGmax_new = max(PGmax_new, PG); 
    562874            PGmin_new = min(PGmin_new, PG); 
    563  
     875             
    564876            if(fabs(PG) > 1.0e-12) 
    565877            { 
     
    575887            } 
    576888        } 
    577  
     889         
    578890        iter++; 
    579891        if(iter % 10 == 0) 
    580         { 
    581892            info("."); 
    582             info_flush(); 
    583         } 
    584  
     893         
    585894        if(PGmax_new - PGmin_new <= eps) 
    586895        { 
     
    590899            { 
    591900                active_size = l; 
    592                 info("*"); info_flush(); 
     901                info("*"); 
    593902                PGmax_old = INF; 
    594903                PGmin_old = -INF; 
     
    603912            PGmin_old = -INF; 
    604913    } 
    605  
     914     
    606915    info("\noptimization finished, #iter = %d\n",iter); 
    607916    if (iter >= max_iter) 
    608         info("Warning: reaching max number of iterations\n"); 
    609  
     917        info("\nWARNING: reaching max number of iterations\nUsing -s 2 may be faster (also see FAQ)\n\n"); 
     918     
    610919    // calculate objective value 
    611  
     920     
    612921    double v = 0; 
    613922    int nSV = 0; 
    614     for(i=0; i<n; i++) 
     923    for(i=0; i<w_size; i++) 
    615924        v += w[i]*w[i]; 
    616925    for(i=0; i<l; i++) 
    617926    { 
    618         if (y[i] == 1) 
    619             v += alpha[i]*(alpha[i]*diag_p - 2); 
    620         else 
    621             v += alpha[i]*(alpha[i]*diag_n - 2); 
     927        v += alpha[i]*(alpha[i]*diag[GETI(i)] - 2); 
    622928        if(alpha[i] > 0) 
    623929            ++nSV; 
     
    625931    info("Objective value = %lf\n",v/2); 
    626932    info("nSV = %d\n",nSV); 
    627  
     933     
    628934    delete [] QD; 
    629935    delete [] alpha; 
     
    632938} 
    633939 
     940// A coordinate descent algorithm for  
     941// the dual of L2-regularized logistic regression problems 
     942// 
     943//  min_\alpha  0.5(\alpha^T Q \alpha) + \sum \alpha_i log (\alpha_i) + (upper_bound_i - alpha_i) log (upper_bound_i - alpha_i) , 
     944//    s.t.      0 <= alpha_i <= upper_bound_i, 
     945//  
     946//  where Qij = yi yj xi^T xj and  
     947//  upper_bound_i = Cp if y_i = 1 
     948//  upper_bound_i = Cn if y_i = -1 
     949// 
     950// Given:  
     951// x, y, Cp, Cn 
     952// eps is the stopping tolerance 
     953// 
     954// solution will be put in w 
     955// 
     956// See Algorithm 5 of Yu et al., MLJ 2010 
     957 
     958#undef GETI 
     959#define GETI(i) (y[i]+1) 
     960// To support weights for instances, use GETI(i) (i) 
     961 
     962void solve_l2r_lr_dual(const problem *prob, double *w, double eps, double Cp, double Cn) 
     963{ 
     964    int l = prob->l; 
     965    int w_size = prob->n; 
     966    int i, s, iter = 0; 
     967    double *xTx = new double[l]; 
     968    int max_iter = 1000; 
     969    int *index = new int[l];         
     970    double *alpha = new double[2*l]; // store alpha and C - alpha 
     971    schar *y = new schar[l];     
     972    int max_inner_iter = 100; // for inner Newton 
     973    double innereps = 1e-2;  
     974    double innereps_min = min(1e-8, eps); 
     975    double upper_bound[3] = {Cn, 0, Cp}; 
     976     
     977    for(i=0; i<w_size; i++) 
     978        w[i] = 0; 
     979    for(i=0; i<l; i++) 
     980    { 
     981        if(prob->y[i] > 0) 
     982        { 
     983            y[i] = +1;  
     984        } 
     985        else 
     986        { 
     987            y[i] = -1; 
     988        } 
     989        alpha[2*i] = min(0.001*upper_bound[GETI(i)], 1e-8); 
     990        alpha[2*i+1] = upper_bound[GETI(i)] - alpha[2*i]; 
     991         
     992        xTx[i] = 0; 
     993        feature_node *xi = prob->x[i]; 
     994        while (xi->index != -1) 
     995        { 
     996            xTx[i] += (xi->value)*(xi->value); 
     997            w[xi->index-1] += y[i]*alpha[2*i]*xi->value; 
     998            xi++; 
     999        } 
     1000        index[i] = i; 
     1001    } 
     1002     
     1003    while (iter < max_iter) 
     1004    { 
     1005        for (i=0; i<l; i++) 
     1006        { 
     1007            int j = i+rand()%(l-i); 
     1008            swap(index[i], index[j]); 
     1009        } 
     1010        int newton_iter = 0; 
     1011        double Gmax = 0; 
     1012        for (s=0; s<l; s++) 
     1013        { 
     1014            i = index[s]; 
     1015            schar yi = y[i]; 
     1016            double C = upper_bound[GETI(i)]; 
     1017            double ywTx = 0, xisq = xTx[i]; 
     1018            feature_node *xi = prob->x[i]; 
     1019            while (xi->index != -1) 
     1020            { 
     1021                ywTx += w[xi->index-1]*xi->value; 
     1022                xi++; 
     1023            } 
     1024            ywTx *= y[i]; 
     1025            double a = xisq, b = ywTx; 
     1026             
     1027            // Decide to minimize g_1(z) or g_2(z) 
     1028            int ind1 = 2*i, ind2 = 2*i+1, sign = 1; 
     1029            if(0.5*a*(alpha[ind2]-alpha[ind1])+b < 0)  
     1030            { 
     1031                ind1 = 2*i+1; 
     1032                ind2 = 2*i; 
     1033                sign = -1; 
     1034            } 
     1035             
     1036            //  g_t(z) = z*log(z) + (C-z)*log(C-z) + 0.5a(z-alpha_old)^2 + sign*b(z-alpha_old) 
     1037            double alpha_old = alpha[ind1]; 
     1038            double z = alpha_old; 
     1039            if(C - z < 0.5 * C)  
     1040                z = 0.1*z; 
     1041            double gp = a*(z-alpha_old)+sign*b+log(z/(C-z)); 
     1042            Gmax = max(Gmax, fabs(gp)); 
     1043             
     1044            // Newton method on the sub-problem 
     1045            const double eta = 0.1; // xi in the paper 
     1046            int inner_iter = 0; 
     1047            while (inner_iter <= max_inner_iter)  
     1048            { 
     1049                if(fabs(gp) < innereps) 
     1050                    break; 
     1051                double gpp = a + C/(C-z)/z; 
     1052                double tmpz = z - gp/gpp; 
     1053                if(tmpz <= 0)  
     1054                    z *= eta; 
     1055                else // tmpz in (0, C) 
     1056                    z = tmpz; 
     1057                gp = a*(z-alpha_old)+sign*b+log(z/(C-z)); 
     1058                newton_iter++; 
     1059                inner_iter++; 
     1060            } 
     1061             
     1062            if(inner_iter > 0) // update w 
     1063            { 
     1064                alpha[ind1] = z; 
     1065                alpha[ind2] = C-z; 
     1066                xi = prob->x[i]; 
     1067                while (xi->index != -1) 
     1068                { 
     1069                    w[xi->index-1] += sign*(z-alpha_old)*yi*xi->value; 
     1070                    xi++; 
     1071                }   
     1072            } 
     1073        } 
     1074         
     1075        iter++; 
     1076        if(iter % 10 == 0) 
     1077            info("."); 
     1078         
     1079        if(Gmax < eps)  
     1080            break; 
     1081         
     1082        if(newton_iter < l/10)  
     1083            innereps = max(innereps_min, 0.1*innereps); 
     1084         
     1085    } 
     1086     
     1087    info("\noptimization finished, #iter = %d\n",iter); 
     1088    if (iter >= max_iter) 
     1089        info("\nWARNING: reaching max number of iterations\nUsing -s 0 may be faster (also see FAQ)\n\n"); 
     1090     
     1091    // calculate objective value 
     1092     
     1093    double v = 0; 
     1094    for(i=0; i<w_size; i++) 
     1095        v += w[i] * w[i]; 
     1096    v *= 0.5; 
     1097    for(i=0; i<l; i++) 
     1098        v += alpha[2*i] * log(alpha[2*i]) + alpha[2*i+1] * log(alpha[2*i+1])  
     1099        - upper_bound[GETI(i)] * log(upper_bound[GETI(i)]); 
     1100    info("Objective value = %lf\n", v); 
     1101     
     1102    delete [] xTx; 
     1103    delete [] alpha; 
     1104    delete [] y; 
     1105    delete [] index; 
     1106} 
     1107 
     1108// A coordinate descent algorithm for  
     1109// L1-regularized L2-loss support vector classification 
     1110// 
     1111//  min_w \sum |wj| + C \sum max(0, 1-yi w^T xi)^2, 
     1112// 
     1113// Given:  
     1114// x, y, Cp, Cn 
     1115// eps is the stopping tolerance 
     1116// 
     1117// solution will be put in w 
     1118// 
     1119// See Yuan et al. (2010) and appendix of LIBLINEAR paper, Fan et al. (2008) 
     1120 
     1121#undef GETI 
     1122#define GETI(i) (y[i]+1) 
     1123// To support weights for instances, use GETI(i) (i) 
     1124 
     1125static void solve_l1r_l2_svc( 
     1126                             problem *prob_col, double *w, double eps,  
     1127                             double Cp, double Cn) 
     1128{ 
     1129    int l = prob_col->l; 
     1130    int w_size = prob_col->n; 
     1131    int j, s, iter = 0; 
     1132    int max_iter = 1000; 
     1133    int active_size = w_size; 
     1134    int max_num_linesearch = 20; 
     1135     
     1136    double sigma = 0.01; 
     1137    double d, G_loss, G, H; 
     1138    double Gmax_old = INF; 
     1139    double Gmax_new; 
     1140    double Gmax_init; 
     1141    double d_old, d_diff; 
     1142    double loss_old, loss_new; 
     1143    double appxcond, cond; 
     1144     
     1145    int *index = new int[w_size]; 
     1146    schar *y = new schar[l]; 
     1147    double *b = new double[l]; // b = 1-ywTx 
     1148    double *xj_sq = new double[w_size]; 
     1149    feature_node *x; 
     1150     
     1151    double C[3] = {Cn,0,Cp}; 
     1152     
     1153    for(j=0; j<l; j++) 
     1154    { 
     1155        b[j] = 1; 
     1156        if(prob_col->y[j] > 0) 
     1157            y[j] = 1; 
     1158        else 
     1159            y[j] = -1; 
     1160    } 
     1161    for(j=0; j<w_size; j++) 
     1162    { 
     1163        w[j] = 0; 
     1164        index[j] = j; 
     1165        xj_sq[j] = 0; 
     1166        x = prob_col->x[j]; 
     1167        while(x->index != -1) 
     1168        { 
     1169            int ind = x->index-1; 
     1170            double val = x->value; 
     1171            x->value *= y[ind]; // x->value stores yi*xij 
     1172            xj_sq[j] += C[GETI(ind)]*val*val; 
     1173            x++; 
     1174        } 
     1175    } 
     1176     
     1177    while(iter < max_iter) 
     1178    { 
     1179        Gmax_new  = 0; 
     1180         
     1181        for(j=0; j<active_size; j++) 
     1182        { 
     1183            int i = j+rand()%(active_size-j); 
     1184            swap(index[i], index[j]); 
     1185        } 
     1186         
     1187        for(s=0; s<active_size; s++) 
     1188        { 
     1189            j = index[s]; 
     1190            G_loss = 0; 
     1191            H = 0; 
     1192             
     1193            x = prob_col->x[j]; 
     1194            while(x->index != -1) 
     1195            { 
     1196                int ind = x->index-1; 
     1197                if(b[ind] > 0) 
     1198                { 
     1199                    double val = x->value; 
     1200                    double tmp = C[GETI(ind)]*val; 
     1201                    G_loss -= tmp*b[ind]; 
     1202                    H += tmp*val; 
     1203                } 
     1204                x++; 
     1205            } 
     1206            G_loss *= 2; 
     1207             
     1208            G = G_loss; 
     1209            H *= 2; 
     1210            H = max(H, 1e-12); 
     1211             
     1212            double Gp = G+1; 
     1213            double Gn = G-1; 
     1214            double violation = 0; 
     1215            if(w[j] == 0) 
     1216            { 
     1217                if(Gp < 0) 
     1218                    violation = -Gp; 
     1219                else if(Gn > 0) 
     1220                    violation = Gn; 
     1221                else if(Gp>Gmax_old/l && Gn<-Gmax_old/l) 
     1222                { 
     1223                    active_size--; 
     1224                    swap(index[s], index[active_size]); 
     1225                    s--; 
     1226                    continue; 
     1227                } 
     1228            } 
     1229            else if(w[j] > 0) 
     1230                violation = fabs(Gp); 
     1231            else 
     1232                violation = fabs(Gn); 
     1233             
     1234            Gmax_new = max(Gmax_new, violation); 
     1235             
     1236            // obtain Newton direction d 
     1237            if(Gp <= H*w[j]) 
     1238                d = -Gp/H; 
     1239            else if(Gn >= H*w[j]) 
     1240                d = -Gn/H; 
     1241            else 
     1242                d = -w[j]; 
     1243             
     1244            if(fabs(d) < 1.0e-12) 
     1245                continue; 
     1246             
     1247            double delta = fabs(w[j]+d)-fabs(w[j]) + G*d; 
     1248            d_old = 0; 
     1249            int num_linesearch; 
     1250            for(num_linesearch=0; num_linesearch < max_num_linesearch; num_linesearch++) 
     1251            { 
     1252                d_diff = d_old - d; 
     1253                cond = fabs(w[j]+d)-fabs(w[j]) - sigma*delta; 
     1254                 
     1255                appxcond = xj_sq[j]*d*d + G_loss*d + cond; 
     1256                if(appxcond <= 0) 
     1257                { 
     1258                    x = prob_col->x[j]; 
     1259                    while(x->index != -1) 
     1260                    { 
     1261                        b[x->index-1] += d_diff*x->value; 
     1262                        x++; 
     1263                    } 
     1264                    break; 
     1265                } 
     1266                 
     1267                if(num_linesearch == 0) 
     1268                { 
     1269                    loss_old = 0; 
     1270                    loss_new = 0; 
     1271                    x = prob_col->x[j]; 
     1272                    while(x->index != -1) 
     1273                    { 
     1274                        int ind = x->index-1; 
     1275                        if(b[ind] > 0) 
     1276                            loss_old += C[GETI(ind)]*b[ind]*b[ind]; 
     1277                        double b_new = b[ind] + d_diff*x->value; 
     1278                        b[ind] = b_new; 
     1279                        if(b_new > 0) 
     1280                            loss_new += C[GETI(ind)]*b_new*b_new; 
     1281                        x++; 
     1282                    } 
     1283                } 
     1284                else 
     1285                { 
     1286                    loss_new = 0; 
     1287                    x = prob_col->x[j]; 
     1288                    while(x->index != -1) 
     1289                    { 
     1290                        int ind = x->index-1; 
     1291                        double b_new = b[ind] + d_diff*x->value; 
     1292                        b[ind] = b_new; 
     1293                        if(b_new > 0) 
     1294                            loss_new += C[GETI(ind)]*b_new*b_new; 
     1295                        x++; 
     1296                    } 
     1297                } 
     1298                 
     1299                cond = cond + loss_new - loss_old; 
     1300                if(cond <= 0) 
     1301                    break; 
     1302                else 
     1303                { 
     1304                    d_old = d; 
     1305                    d *= 0.5; 
     1306                    delta *= 0.5; 
     1307                } 
     1308            } 
     1309             
     1310            w[j] += d; 
     1311             
     1312            // recompute b[] if line search takes too many steps 
     1313            if(num_linesearch >= max_num_linesearch) 
     1314            { 
     1315                info("#"); 
     1316                for(int i=0; i<l; i++) 
     1317                    b[i] = 1; 
     1318                 
     1319                for(int i=0; i<w_size; i++) 
     1320                { 
     1321                    if(w[i]==0) continue; 
     1322                    x = prob_col->x[i]; 
     1323                    while(x->index != -1) 
     1324                    { 
     1325                        b[x->index-1] -= w[i]*x->value; 
     1326                        x++; 
     1327                    } 
     1328                } 
     1329            } 
     1330        } 
     1331         
     1332        if(iter == 0) 
     1333            Gmax_init = Gmax_new; 
     1334        iter++; 
     1335        if(iter % 10 == 0) 
     1336            info("."); 
     1337         
     1338        if(Gmax_new <= eps*Gmax_init) 
     1339        { 
     1340            if(active_size == w_size) 
     1341                break; 
     1342            else 
     1343            { 
     1344                active_size = w_size; 
     1345                info("*"); 
     1346                Gmax_old = INF; 
     1347                continue; 
     1348            } 
     1349        } 
     1350         
     1351        Gmax_old = Gmax_new; 
     1352    } 
     1353     
     1354    info("\noptimization finished, #iter = %d\n", iter); 
     1355    if(iter >= max_iter) 
     1356        info("\nWARNING: reaching max number of iterations\n"); 
     1357     
     1358    // calculate objective value 
     1359     
     1360    double v = 0; 
     1361    int nnz = 0; 
     1362    for(j=0; j<w_size; j++) 
     1363    { 
     1364        x = prob_col->x[j]; 
     1365        while(x->index != -1) 
     1366        { 
     1367            x->value *= prob_col->y[x->index-1]; // restore x->value 
     1368            x++; 
     1369        } 
     1370        if(w[j] != 0) 
     1371        { 
     1372            v += fabs(w[j]); 
     1373            nnz++; 
     1374        } 
     1375    } 
     1376    for(j=0; j<l; j++) 
     1377        if(b[j] > 0) 
     1378            v += C[GETI(j)]*b[j]*b[j]; 
     1379     
     1380    info("Objective value = %lf\n", v); 
     1381    info("#nonzeros/#features = %d/%d\n", nnz, w_size); 
     1382     
     1383    delete [] index; 
     1384    delete [] y; 
     1385    delete [] b; 
     1386    delete [] xj_sq; 
     1387} 
     1388 
     1389// A coordinate descent algorithm for  
     1390// L1-regularized logistic regression problems 
     1391// 
     1392//  min_w \sum |wj| + C \sum log(1+exp(-yi w^T xi)), 
     1393// 
     1394// Given:  
     1395// x, y, Cp, Cn 
     1396// eps is the stopping tolerance 
     1397// 
     1398// solution will be put in w 
     1399// 
     1400// See Yuan et al. (2010) and appendix of LIBLINEAR paper, Fan et al. (2008) 
     1401 
     1402#undef GETI 
     1403#define GETI(i) (y[i]+1) 
     1404// To support weights for instances, use GETI(i) (i) 
     1405 
     1406static void solve_l1r_lr( 
     1407                         const problem *prob_col, double *w, double eps,  
     1408                         double Cp, double Cn) 
     1409{ 
     1410    int l = prob_col->l; 
     1411    int w_size = prob_col->n; 
     1412    int j, s, iter = 0; 
     1413    int max_iter = 1000; 
     1414    int active_size = w_size; 
     1415    int max_num_linesearch = 20; 
     1416     
     1417    double x_min = 0; 
     1418    double sigma = 0.01; 
     1419    double d, G, H; 
     1420    double Gmax_old = INF; 
     1421    double Gmax_new; 
     1422    double Gmax_init; 
     1423    double sum1, appxcond1; 
     1424    double sum2, appxcond2; 
     1425    double cond; 
     1426     
     1427    int *index = new int[w_size]; 
     1428    schar *y = new schar[l]; 
     1429    double *exp_wTx = new double[l]; 
     1430    double *exp_wTx_new = new double[l]; 
     1431    double *xj_max = new double[w_size]; 
     1432    double *C_sum = new double[w_size]; 
     1433    double *xjneg_sum = new double[w_size]; 
     1434    double *xjpos_sum = new double[w_size]; 
     1435    feature_node *x; 
     1436     
     1437    double C[3] = {Cn,0,Cp}; 
     1438     
     1439    for(j=0; j<l; j++) 
     1440    { 
     1441        exp_wTx[j] = 1; 
     1442        if(prob_col->y[j] > 0) 
     1443            y[j] = 1; 
     1444        else 
     1445            y[j] = -1; 
     1446    } 
     1447    for(j=0; j<w_size; j++) 
     1448    { 
     1449        w[j] = 0; 
     1450        index[j] = j; 
     1451        xj_max[j] = 0; 
     1452        C_sum[j] = 0; 
     1453        xjneg_sum[j] = 0; 
     1454        xjpos_sum[j] = 0; 
     1455        x = prob_col->x[j]; 
     1456        while(x->index != -1) 
     1457        { 
     1458            int ind = x->index-1; 
     1459            double val = x->value; 
     1460            x_min = min(x_min, val); 
     1461            xj_max[j] = max(xj_max[j], val); 
     1462            C_sum[j] += C[GETI(ind)]; 
     1463            if(y[ind] == -1) 
     1464                xjneg_sum[j] += C[GETI(ind)]*val; 
     1465            else 
     1466                xjpos_sum[j] += C[GETI(ind)]*val; 
     1467            x++; 
     1468        } 
     1469    } 
     1470     
     1471    while(iter < max_iter) 
     1472    { 
     1473        Gmax_new = 0; 
     1474         
     1475        for(j=0; j<active_size; j++) 
     1476        { 
     1477            int i = j+rand()%(active_size-j); 
     1478            swap(index[i], index[j]); 
     1479        } 
     1480         
     1481        for(s=0; s<active_size; s++) 
     1482        { 
     1483            j = index[s]; 
     1484            sum1 = 0; 
     1485            sum2 = 0; 
     1486            H = 0; 
     1487             
     1488            x = prob_col->x[j]; 
     1489            while(x->index != -1) 
     1490            { 
     1491                int ind = x->index-1; 
     1492                double exp_wTxind = exp_wTx[ind]; 
     1493                double tmp1 = x->value/(1+exp_wTxind); 
     1494                double tmp2 = C[GETI(ind)]*tmp1; 
     1495                double tmp3 = tmp2*exp_wTxind; 
     1496                sum2 += tmp2; 
     1497                sum1 += tmp3; 
     1498                H += tmp1*tmp3; 
     1499                x++; 
     1500            } 
     1501             
     1502            G = -sum2 + xjneg_sum[j]; 
     1503             
     1504            double Gp = G+1; 
     1505            double Gn = G-1; 
     1506            double violation = 0; 
     1507            if(w[j] == 0) 
     1508            { 
     1509                if(Gp < 0) 
     1510                    violation = -Gp; 
     1511                else if(Gn > 0) 
     1512                    violation = Gn; 
     1513                else if(Gp>Gmax_old/l && Gn<-Gmax_old/l) 
     1514                { 
     1515                    active_size--; 
     1516                    swap(index[s], index[active_size]); 
     1517                    s--; 
     1518                    continue; 
     1519                } 
     1520            } 
     1521            else if(w[j] > 0) 
     1522                violation = fabs(Gp); 
     1523            else 
     1524                violation = fabs(Gn); 
     1525             
     1526            Gmax_new = max(Gmax_new, violation); 
     1527             
     1528            // obtain Newton direction d 
     1529            if(Gp <= H*w[j]) 
     1530                d = -Gp/H; 
     1531            else if(Gn >= H*w[j]) 
     1532                d = -Gn/H; 
     1533            else 
     1534                d = -w[j]; 
     1535             
     1536            if(fabs(d) < 1.0e-12) 
     1537                continue; 
     1538             
     1539            d = min(max(d,-10.0),10.0); 
     1540             
     1541            double delta = fabs(w[j]+d)-fabs(w[j]) + G*d; 
     1542            int num_linesearch; 
     1543            for(num_linesearch=0; num_linesearch < max_num_linesearch; num_linesearch++) 
     1544            { 
     1545                cond = fabs(w[j]+d)-fabs(w[j]) - sigma*delta; 
     1546                 
     1547                if(x_min >= 0) 
     1548                { 
     1549                    double tmp = exp(d*xj_max[j]); 
     1550                    appxcond1 = log(1+sum1*(tmp-1)/xj_max[j]/C_sum[j])*C_sum[j] + cond - d*xjpos_sum[j]; 
     1551                    appxcond2 = log(1+sum2*(1/tmp-1)/xj_max[j]/C_sum[j])*C_sum[j] + cond + d*xjneg_sum[j]; 
     1552                    if(min(appxcond1,appxcond2) <= 0) 
     1553                    { 
     1554                        x = prob_col->x[j]; 
     1555                        while(x->index != -1) 
     1556                        { 
     1557                            exp_wTx[x->index-1] *= exp(d*x->value); 
     1558                            x++; 
     1559                        } 
     1560                        break; 
     1561                    } 
     1562                } 
     1563                 
     1564                cond += d*xjneg_sum[j]; 
     1565                 
     1566                int i = 0; 
     1567                x = prob_col->x[j]; 
     1568                while(x->index != -1) 
     1569                { 
     1570                    int ind = x->index-1; 
     1571                    double exp_dx = exp(d*x->value); 
     1572                    exp_wTx_new[i] = exp_wTx[ind]*exp_dx; 
     1573                    cond += C[GETI(ind)]*log((1+exp_wTx_new[i])/(exp_dx+exp_wTx_new[i])); 
     1574                    x++; i++; 
     1575                } 
     1576                 
     1577                if(cond <= 0) 
     1578                { 
     1579                    int i = 0; 
     1580                    x = prob_col->x[j]; 
     1581                    while(x->index != -1) 
     1582                    { 
     1583                        int ind = x->index-1; 
     1584                        exp_wTx[ind] = exp_wTx_new[i]; 
     1585                        x++; i++; 
     1586                    } 
     1587                    break; 
     1588                } 
     1589                else 
     1590                { 
     1591                    d *= 0.5; 
     1592                    delta *= 0.5; 
     1593                } 
     1594            } 
     1595             
     1596            w[j] += d; 
     1597             
     1598            // recompute exp_wTx[] if line search takes too many steps 
     1599            if(num_linesearch >= max_num_linesearch) 
     1600            { 
     1601                info("#"); 
     1602                for(int i=0; i<l; i++) 
     1603                    exp_wTx[i] = 0; 
     1604                 
     1605                for(int i=0; i<w_size; i++) 
     1606                { 
     1607                    if(w[i]==0) continue; 
     1608                    x = prob_col->x[i]; 
     1609                    while(x->index != -1) 
     1610                    { 
     1611                        exp_wTx[x->index-1] += w[i]*x->value; 
     1612                        x++; 
     1613                    } 
     1614                } 
     1615                 
     1616                for(int i=0; i<l; i++) 
     1617                    exp_wTx[i] = exp(exp_wTx[i]); 
     1618            } 
     1619        } 
     1620         
     1621        if(iter == 0) 
     1622            Gmax_init = Gmax_new; 
     1623        iter++; 
     1624        if(iter % 10 == 0) 
     1625            info("."); 
     1626         
     1627        if(Gmax_new <= eps*Gmax_init) 
     1628        { 
     1629            if(active_size == w_size) 
     1630                break; 
     1631            else 
     1632            { 
     1633                active_size = w_size; 
     1634                info("*"); 
     1635                Gmax_old = INF; 
     1636                continue; 
     1637            } 
     1638        } 
     1639         
     1640        Gmax_old = Gmax_new; 
     1641    } 
     1642     
     1643    info("\noptimization finished, #iter = %d\n", iter); 
     1644    if(iter >= max_iter) 
     1645        info("\nWARNING: reaching max number of iterations\n"); 
     1646     
     1647    // calculate objective value 
     1648     
     1649    double v = 0; 
     1650    int nnz = 0; 
     1651    for(j=0; j<w_size; j++) 
     1652        if(w[j] != 0) 
     1653        { 
     1654            v += fabs(w[j]); 
     1655            nnz++; 
     1656        } 
     1657    for(j=0; j<l; j++) 
     1658        if(y[j] == 1) 
     1659            v += C[GETI(j)]*log(1+1/exp_wTx[j]); 
     1660        else 
     1661            v += C[GETI(j)]*log(1+exp_wTx[j]); 
     1662     
     1663    info("Objective value = %lf\n", v); 
     1664    info("#nonzeros/#features = %d/%d\n", nnz, w_size); 
     1665     
     1666    delete [] index; 
     1667    delete [] y; 
     1668    delete [] exp_wTx; 
     1669    delete [] exp_wTx_new; 
     1670    delete [] xj_max; 
     1671    delete [] C_sum; 
     1672    delete [] xjneg_sum; 
     1673    delete [] xjpos_sum; 
     1674} 
     1675 
     1676// transpose matrix X from row format to column format 
     1677static void transpose(const problem *prob, feature_node **x_space_ret, problem *prob_col) 
     1678{ 
     1679    int i; 
     1680    int l = prob->l; 
     1681    int n = prob->n; 
     1682    int nnz = 0; 
     1683    int *col_ptr = new int[n+1]; 
     1684    feature_node *x_space; 
     1685    prob_col->l = l; 
     1686    prob_col->n = n; 
     1687    prob_col->y = new int[l]; 
     1688    prob_col->x = new feature_node*[n]; 
     1689     
     1690    for(i=0; i<l; i++) 
     1691        prob_col->y[i] = prob->y[i]; 
     1692     
     1693    for(i=0; i<n+1; i++) 
     1694        col_ptr[i] = 0; 
     1695    for(i=0; i<l; i++) 
     1696    { 
     1697        feature_node *x = prob->x[i]; 
     1698        while(x->index != -1) 
     1699        { 
     1700            nnz++; 
     1701            col_ptr[x->index]++; 
     1702            x++; 
     1703        } 
     1704    } 
     1705    for(i=1; i<n+1; i++) 
     1706        col_ptr[i] += col_ptr[i-1] + 1; 
     1707     
     1708    x_space = new feature_node[nnz+n]; 
     1709    for(i=0; i<n; i++) 
     1710        prob_col->x[i] = &x_space[col_ptr[i]]; 
     1711     
     1712    for(i=0; i<l; i++) 
     1713    { 
     1714        feature_node *x = prob->x[i]; 
     1715        while(x->index != -1) 
     1716        { 
     1717            int ind = x->index-1; 
     1718            x_space[col_ptr[ind]].index = i+1; // starts from 1 
     1719            x_space[col_ptr[ind]].value = x->value; 
     1720            col_ptr[ind]++; 
     1721            x++; 
     1722        } 
     1723    } 
     1724    for(i=0; i<n; i++) 
     1725        x_space[col_ptr[i]].index = -1; 
     1726     
     1727    *x_space_ret = x_space; 
     1728     
     1729    delete [] col_ptr; 
     1730} 
     1731 
    6341732// label: label name, start: begin of each class, count: #data of classes, perm: indices to the original data 
    6351733// perm, length l, must be allocated before calling this subroutine 
    636 void group_classes(const problem *prob, int *nr_class_ret, int **label_ret, int **start_ret, int **count_ret, int *perm) 
     1734static void group_classes(const problem *prob, int *nr_class_ret, int **label_ret, int **start_ret, int **count_ret, int *perm) 
    6371735{ 
    6381736    int l = prob->l; 
     
    6431741    int *data_label = Malloc(int,l); 
    6441742    int i; 
    645  
     1743     
    6461744    for(i=0;i<l;i++) 
    6471745    { 
     
    6701768        } 
    6711769    } 
    672  
     1770     
    6731771    int *start = Malloc(int,nr_class); 
    6741772    start[0] = 0; 
     
    6831781    for(i=1;i<nr_class;i++) 
    6841782        start[i] = start[i-1]+count[i-1]; 
    685  
     1783     
    6861784    *nr_class_ret = nr_class; 
    6871785    *label_ret = label; 
     
    6911789} 
    6921790 
    693 void train_one(const problem *prob, const parameter *param, double *w, double Cp, double Cn) 
     1791static void train_one(const problem *prob, const parameter *param, double *w, double Cp, double Cn) 
    6941792{ 
    6951793    double eps=param->eps; 
    6961794    int pos = 0; 
    6971795    int neg = 0; 
    698     for (int i=0; i<prob->l;i++) 
    699         if (prob->y[i]==+1) 
     1796    for(int i=0;i<prob->l;i++) 
     1797        if(prob->y[i]==+1) 
    7001798            pos++; 
    7011799    neg = prob->l - pos; 
    702  
    703     function1 *fun_obj=NULL; 
     1800     
     1801    function *fun_obj=NULL; 
    7041802    switch(param->solver_type) 
    7051803    { 
    706         case L2_LR: 
    707         { 
    708             fun_obj=new l2_lr_fun(prob, Cp, Cn); 
     1804        case L2R_LR: 
     1805        { 
     1806            fun_obj=new l2r_lr_fun(prob, Cp, Cn); 
    7091807            TRON tron_obj(fun_obj, eps*min(pos,neg)/prob->l); 
     1808            tron_obj.set_print_string(liblinear_print_string); 
    7101809            tron_obj.tron(w); 
    7111810            delete fun_obj; 
    7121811            break; 
    7131812        } 
    714         case L2LOSS_SVM: 
    715         { 
    716             fun_obj=new l2loss_svm_fun(prob, Cp, Cn); 
     1813        case L2R_L2LOSS_SVC: 
     1814        { 
     1815            fun_obj=new l2r_l2_svc_fun(prob, Cp, Cn); 
    7171816            TRON tron_obj(fun_obj, eps*min(pos,neg)/prob->l); 
     1817            tron_obj.set_print_string(liblinear_print_string); 
    7181818            tron_obj.tron(w); 
    7191819            delete fun_obj; 
    7201820            break; 
    7211821        } 
    722         case L2LOSS_SVM_DUAL: 
    723             solve_linear_c_svc(prob, w, eps, Cp, Cn, L2LOSS_SVM_DUAL); 
     1822        case L2R_L2LOSS_SVC_DUAL: 
     1823            solve_l2r_l1l2_svc(prob, w, eps, Cp, Cn, L2R_L2LOSS_SVC_DUAL); 
    7241824            break; 
    725         case L1LOSS_SVM_DUAL: 
    726             solve_linear_c_svc(prob, w, eps, Cp, Cn, L1LOSS_SVM_DUAL); 
     1825        case L2R_L1LOSS_SVC_DUAL: 
     1826            solve_l2r_l1l2_svc(prob, w, eps, Cp, Cn, L2R_L1LOSS_SVC_DUAL); 
     1827            break; 
     1828        case L1R_L2LOSS_SVC: 
     1829        { 
     1830            problem prob_col; 
     1831            feature_node *x_space = NULL; 
     1832            transpose(prob, &x_space ,&prob_col); 
     1833            solve_l1r_l2_svc(&prob_col, w, eps*min(pos,neg)/prob->l, Cp, Cn); 
     1834            delete [] prob_col.y; 
     1835            delete [] prob_col.x; 
     1836            delete [] x_space; 
     1837            break; 
     1838        } 
     1839        case L1R_LR: 
     1840        { 
     1841            problem prob_col; 
     1842            feature_node *x_space = NULL; 
     1843            transpose(prob, &x_space ,&prob_col); 
     1844            solve_l1r_lr(&prob_col, w, eps*min(pos,neg)/prob->l, Cp, Cn); 
     1845            delete [] prob_col.y; 
     1846            delete [] prob_col.x; 
     1847            delete [] x_space; 
     1848            break; 
     1849        } 
     1850        case L2R_LR_DUAL: 
     1851            solve_l2r_lr_dual(prob, w, eps, Cp, Cn); 
    7271852            break; 
    7281853        default: 
     
    7331858 
    7341859// 
    735 // Interface function1s 
     1860// Interface functions 
    7361861// 
    7371862model* train(const problem *prob, const parameter *param) 
    7381863{ 
    739     int i; 
     1864    int i,j; 
    7401865    int l = prob->l; 
    7411866    int n = prob->n; 
     1867    int w_size = prob->n; 
    7421868    model *model_ = Malloc(model,1); 
    743  
     1869     
    7441870    if(prob->bias>=0) 
    7451871        model_->nr_feature=n-1; 
     
    7481874    model_->param = *param; 
    7491875    model_->bias = prob->bias; 
    750  
     1876     
    7511877    int nr_class; 
    7521878    int *label = NULL; 
     
    7541880    int *count = NULL; 
    7551881    int *perm = Malloc(int,l); 
    756  
     1882     
    7571883    // group training data of the same class 
    7581884    group_classes(prob,&nr_class,&label,&start,&count,perm); 
    759  
     1885     
    7601886    model_->nr_class=nr_class; 
    7611887    model_->label = Malloc(int,nr_class); 
    7621888    for(i=0;i<nr_class;i++) 
    7631889        model_->label[i] = label[i]; 
    764  
     1890     
    7651891    // calculate weighted C 
    7661892    double *weighted_C = Malloc(double, nr_class); 
     
    7691895    for(i=0;i<param->nr_weight;i++) 
    7701896    { 
    771         int j; 
    7721897        for(j=0;j<nr_class;j++) 
    7731898            if(param->weight_label[i] == label[j]) 
    7741899                break; 
    7751900        if(j == nr_class) 
    776             fprintf(stderr,"warning: class label %d specified in weight is not found\n", param->weight_label[i]); 
     1901            fprintf(stderr,"WARNING: class label %d specified in weight is not found\n", param->weight_label[i]); 
    7771902        else 
    7781903            weighted_C[j] *= param->weight[i]; 
    7791904    } 
    780  
     1905     
    7811906    // constructing the subproblem 
    7821907    feature_node **x = Malloc(feature_node *,l); 
    7831908    for(i=0;i<l;i++) 
    7841909        x[i] = prob->x[perm[i]]; 
    785  
     1910     
    7861911    int k; 
    7871912    problem sub_prob; 
     
    7901915    sub_prob.x = Malloc(feature_node *,sub_prob.l); 
    7911916    sub_prob.y = Malloc(int,sub_prob.l); 
    792  
     1917     
    7931918    for(k=0; k<sub_prob.l; k++) 
    7941919        sub_prob.x[k] = x[k]; 
    795  
    796     if(nr_class==2) 
    797     { 
    798         model_->w=Malloc(double, n); 
    799  
    800         int e0 = start[0]+count[0]; 
    801         k=0; 
    802         for(; k<e0; k++) 
    803             sub_prob.y[k] = +1; 
    804         for(; k<sub_prob.l; k++) 
    805             sub_prob.y[k] = -1; 
    806  
    807         train_one(&sub_prob, param, &model_->w[0], weighted_C[0], weighted_C[1]); 
     1920     
     1921    // multi-class svm by Crammer and Singer 
     1922    if(param->solver_type == MCSVM_CS) 
     1923    { 
     1924        model_->w=Malloc(double, n*nr_class); 
     1925        for(i=0;i<nr_class;i++) 
     1926            for(j=start[i];j<start[i]+count[i];j++) 
     1927                sub_prob.y[j] = i; 
     1928        Solver_MCSVM_CS Solver(&sub_prob, nr_class, weighted_C, param->eps); 
     1929        Solver.Solve(model_->w); 
    8081930    } 
    8091931    else 
    8101932    { 
    811         model_->w=Malloc(double, n*nr_class); 
    812         double *w=Malloc(double, n); 
    813         for(i=0;i<nr_class;i++) 
    814         { 
    815             int si = start[i]; 
    816             int ei = si+count[i]; 
    817  
     1933        if(nr_class == 2) 
     1934        { 
     1935            model_->w=Malloc(double, w_size); 
     1936             
     1937            int e0 = start[0]+count[0]; 
    8181938            k=0; 
    819             for(; k<si; k++) 
    820                 sub_prob.y[k] = -1; 
    821             for(; k<ei; k++) 
     1939            for(; k<e0; k++) 
    8221940                sub_prob.y[k] = +1; 
    8231941            for(; k<sub_prob.l; k++) 
    8241942                sub_prob.y[k] = -1; 
    825  
    826             train_one(&sub_prob, param, w, weighted_C[i], param->C); 
    827  
    828             for(int j=0;j<n;j++) 
    829                 model_->w[j*nr_class+i] = w[j]; 
    830         } 
    831         free(w); 
    832     } 
    833  
     1943             
     1944            train_one(&sub_prob, param, &model_->w[0], weighted_C[0], weighted_C[1]); 
     1945        } 
     1946        else 
     1947        { 
     1948            model_->w=Malloc(double, w_size*nr_class); 
     1949            double *w=Malloc(double, w_size); 
     1950            for(i=0;i<nr_class;i++) 
     1951            { 
     1952                int si = start[i]; 
     1953                int ei = si+count[i]; 
     1954                 
     1955                k=0; 
     1956                for(; k<si; k++) 
     1957                    sub_prob.y[k] = -1; 
     1958                for(; k<ei; k++) 
     1959                    sub_prob.y[k] = +1; 
     1960                for(; k<sub_prob.l; k++) 
     1961                    sub_prob.y[k] = -1; 
     1962                 
     1963                train_one(&sub_prob, param, w, weighted_C[i], param->C); 
     1964                 
     1965                for(int j=0;j<w_size;j++) 
     1966                    model_->w[j*nr_class+i] = w[j]; 
     1967            } 
     1968            free(w); 
     1969        } 
     1970         
     1971    } 
     1972     
    8341973    free(x); 
    8351974    free(label); 
     
    8431982} 
    8441983 
    845 void destroy_model(struct model *model_) 
    846 { 
    847     if(model_->w != NULL) 
    848         free(model_->w); 
    849     if(model_->label != NULL) 
    850         free(model_->label); 
    851     free(model_); 
    852 } 
    853  
    854 const char *solver_type_table[]= 
    855 { 
    856     "L2_LR", "L2LOSS_SVM_DUAL", "L2LOSS_SVM","L1LOSS_SVM_DUAL", NULL 
     1984void cross_validation(const problem *prob, const parameter *param, int nr_fold, int *target) 
     1985{ 
     1986    int i; 
     1987    int *fold_start = Malloc(int,nr_fold+1); 
     1988    int l = prob->l; 
     1989    int *perm = Malloc(int,l); 
     1990     
     1991    for(i=0;i<l;i++) perm[i]=i; 
     1992    for(i=0;i<l;i++) 
     1993    { 
     1994        int j = i+rand()%(l-i); 
     1995        swap(perm[i],perm[j]); 
     1996    } 
     1997    for(i=0;i<=nr_fold;i++) 
     1998        fold_start[i]=i*l/nr_fold; 
     1999     
     2000    for(i=0;i<nr_fold;i++) 
     2001    { 
     2002        int begin = fold_start[i]; 
     2003        int end = fold_start[i+1]; 
     2004        int j,k; 
     2005        struct problem subprob; 
     2006         
     2007        subprob.bias = prob->bias; 
     2008        subprob.n = prob->n; 
     2009        subprob.l = l-(end-begin); 
     2010        subprob.x = Malloc(struct feature_node*,subprob.l); 
     2011        subprob.y = Malloc(int,subprob.l); 
     2012         
     2013        k=0; 
     2014        for(j=0;j<begin;j++) 
     2015        { 
     2016            subprob.x[k] = prob->x[perm[j]]; 
     2017            subprob.y[k] = prob->y[perm[j]]; 
     2018            ++k; 
     2019        } 
     2020        for(j=end;j<l;j++) 
     2021        { 
     2022            subprob.x[k] = prob->x[perm[j]]; 
     2023            subprob.y[k] = prob->y[perm[j]]; 
     2024            ++k; 
     2025        } 
     2026        struct model *submodel = train(&subprob,param); 
     2027        for(j=begin;j<end;j++) 
     2028            target[perm[j]] = predict(submodel,prob->x[perm[j]]); 
     2029        free_and_destroy_model(&submodel); 
     2030        free(subprob.x); 
     2031        free(subprob.y); 
     2032    } 
     2033    free(fold_start); 
     2034    free(perm); 
     2035} 
     2036 
     2037int predict_values(const struct model *model_, const struct feature_node *x, double *dec_values) 
     2038{ 
     2039    int idx; 
     2040    int n; 
     2041    if(model_->bias>=0) 
     2042        n=model_->nr_feature+1; 
     2043    else 
     2044        n=model_->nr_feature; 
     2045    double *w=model_->w; 
     2046    int nr_class=model_->nr_class; 
     2047    int i; 
     2048    int nr_w; 
     2049    if(nr_class==2 && model_->param.solver_type != MCSVM_CS) 
     2050        nr_w = 1; 
     2051    else 
     2052        nr_w = nr_class; 
     2053     
     2054    const feature_node *lx=x; 
     2055    for(i=0;i<nr_w;i++) 
     2056        dec_values[i] = 0; 
     2057    for(; (idx=lx->index)!=-1; lx++) 
     2058    { 
     2059        // the dimension of testing data may exceed that of training 
     2060        if(idx<=n) 
     2061            for(i=0;i<nr_w;i++) 
     2062                dec_values[i] += w[(idx-1)*nr_w+i]*lx->value; 
     2063    } 
     2064     
     2065    if(nr_class==2) 
     2066        return (dec_values[0]>0)?model_->label[0]:model_->label[1]; 
     2067    else 
     2068    { 
     2069        int dec_max_idx = 0; 
     2070        for(i=1;i<nr_class;i++) 
     2071        { 
     2072            if(dec_values[i] > dec_values[dec_max_idx]) 
     2073                dec_max_idx = i; 
     2074        } 
     2075        return model_->label[dec_max_idx]; 
     2076    } 
     2077} 
     2078 
     2079int predict(const model *model_, const feature_node *x) 
     2080{ 
     2081    double *dec_values = Malloc(double, model_->nr_class); 
     2082    int label=predict_values(model_, x, dec_values); 
     2083    free(dec_values); 
     2084    return label; 
     2085} 
     2086 
     2087int predict_probability(const struct model *model_, const struct feature_node *x, double* prob_estimates) 
     2088{ 
     2089    if(check_probability_model(model_)) 
     2090    { 
     2091        int i; 
     2092        int nr_class=model_->nr_class; 
     2093        int nr_w; 
     2094        if(nr_class==2) 
     2095            nr_w = 1; 
     2096        else 
     2097            nr_w = nr_class; 
     2098         
     2099        int label=predict_values(model_, x, prob_estimates); 
     2100        for(i=0;i<nr_w;i++) 
     2101            prob_estimates[i]=1/(1+exp(-prob_estimates[i])); 
     2102         
     2103        if(nr_class==2) // for binary classification 
     2104            prob_estimates[1]=1.-prob_estimates[0]; 
     2105        else 
     2106        { 
     2107            double sum=0; 
     2108            for(i=0; i<nr_class; i++) 
     2109                sum+=prob_estimates[i]; 
     2110             
     2111            for(i=0; i<nr_class; i++) 
     2112                prob_estimates[i]=prob_estimates[i]/sum; 
     2113        } 
     2114         
     2115        return label;        
     2116    } 
     2117    else 
     2118        return 0; 
     2119} 
     2120 
     2121static const char *solver_type_table[]= 
     2122{ 
     2123    "L2R_LR", "L2R_L2LOSS_SVC_DUAL", "L2R_L2LOSS_SVC", "L2R_L1LOSS_SVC_DUAL", "MCSVM_CS", 
     2124    "L1R_L2LOSS_SVC", "L1R_LR", "L2R_LR_DUAL", NULL 
    8572125}; 
    8582126 
     
    8632131    int n; 
    8642132    const parameter& param = model_->param; 
    865  
     2133     
    8662134    if(model_->bias>=0) 
    8672135        n=nr_feature+1; 
    8682136    else 
    8692137        n=nr_feature; 
     2138    int w_size = n; 
    8702139    FILE *fp = fopen(model_file_name,"w"); 
    8712140    if(fp==NULL) return -1; 
    872  
    873     int nr_classifier; 
    874     if(model_->nr_class==2) 
    875         nr_classifier=1; 
     2141     
     2142    int nr_w; 
     2143    if(model_->nr_class==2 && model_->param.solver_type != MCSVM_CS) 
     2144        nr_w=1; 
    8762145    else 
    877         nr_classifier=model_->nr_class; 
    878  
     2146        nr_w=model_->nr_class; 
     2147     
    8792148    fprintf(fp, "solver_type %s\n", solver_type_table[param.solver_type]); 
    8802149    fprintf(fp, "nr_class %d\n", model_->nr_class); 
     
    8832152        fprintf(fp, " %d", model_->label[i]); 
    8842153    fprintf(fp, "\n"); 
    885  
     2154     
    8862155    fprintf(fp, "nr_feature %d\n", nr_feature); 
    887  
     2156     
    8882157    fprintf(fp, "bias %.16g\n", model_->bias); 
    889  
     2158     
    8902159    fprintf(fp, "w\n"); 
    891     for(i=0; i<n; i++) 
     2160    for(i=0; i<w_size; i++) 
    8922161    { 
    8932162        int j; 
    894         for(j=0; j<nr_classifier; j++) 
    895             fprintf(fp, "%.16g ", model_->w[i*nr_classifier+j]); 
     2163        for(j=0; j<nr_w; j++) 
     2164            fprintf(fp, "%.16g ", model_->w[i*nr_w+j]); 
    8962165        fprintf(fp, "\n"); 
    8972166    } 
    898  
     2167     
    8992168    if (ferror(fp) != 0 || fclose(fp) != 0) return -1; 
    9002169    else return 0; 
     
    9052174    FILE *fp = fopen(model_file_name,"r"); 
    9062175    if(fp==NULL) return NULL; 
    907  
     2176     
    9082177    int i; 
    9092178    int nr_feature; 
     
    9132182    model *model_ = Malloc(model,1); 
    9142183    parameter& param = model_->param; 
    915  
     2184     
    9162185    model_->label = NULL; 
    917  
     2186     
    9182187    char cmd[81]; 
    9192188    while(1) 
     
    9732242        } 
    9742243    } 
    975  
     2244     
    9762245    nr_feature=model_->nr_feature; 
    9772246    if(model_->bias>=0) 
     
    9792248    else 
    9802249        n=nr_feature; 
    981  
    982     int nr_classifier; 
    983     if(nr_class==2) 
    984         nr_classifier = 1; 
     2250    int w_size = n; 
     2251    int nr_w; 
     2252    if(nr_class==2 && param.solver_type != MCSVM_CS) 
     2253        nr_w = 1; 
    9852254    else 
    986         nr_classifier = nr_class; 
    987  
    988     model_->w=Malloc(double, n*nr_classifier); 
    989     for(i=0; i<n; i++) 
     2255        nr_w = nr_class; 
     2256     
     2257    model_->w=Malloc(double, w_size*nr_w); 
     2258    for(i=0; i<w_size; i++) 
    9902259    { 
    9912260        int j; 
    992         for(j=0; j<nr_classifier; j++) 
    993             fscanf(fp, "%lf ", &model_->w[i*nr_classifier+j]); 
     2261        for(j=0; j<nr_w; j++) 
     2262            fscanf(fp, "%lf ", &model_->w[i*nr_w+j]); 
    9942263        fscanf(fp, "\n"); 
    9952264    } 
    9962265    if (ferror(fp) != 0 || fclose(fp) != 0) return NULL; 
    997  
     2266     
    9982267    return model_; 
    9992268} 
    10002269 
    1001 int predict_values(const struct model *model_, const struct feature_node *x, double *dec_values) 
    1002 { 
    1003     int idx; 
    1004     int n; 
    1005     if(model_->bias>=0) 
    1006         n=model_->nr_feature+1; 
    1007     else 
    1008         n=model_->nr_feature; 
    1009     double *w=model_->w; 
    1010     int nr_class=model_->nr_class; 
    1011     int i; 
    1012     int nr_classifier; 
    1013     if(nr_class==2) 
    1014         nr_classifier = 1; 
    1015     else 
    1016         nr_classifier = nr_class; 
    1017  
    1018     const feature_node *lx=x; 
    1019     for(i=0;i<nr_classifier;i++) 
    1020         dec_values[i] = 0; 
    1021     for(; (idx=lx->index)!=-1; lx++) 
    1022     { 
    1023         // the dimension of testing data may exceed that of training 
    1024         if(idx<=n) 
    1025             for(i=0;i<nr_classifier;i++) 
    1026                 dec_values[i] += w[(idx-1)*nr_classifier+i]*lx->value; 
    1027     } 
    1028  
    1029     if(nr_class==2) 
    1030         return (dec_values[0]>0)?model_->label[0]:model_->label[1]; 
    1031     else 
    1032     { 
    1033         int dec_max_idx = 0; 
    1034         for(i=1;i<nr_class;i++) 
    1035         { 
    1036             if(dec_values[i] > dec_values[dec_max_idx]) 
    1037                 dec_max_idx = i; 
    1038         } 
    1039         return model_->label[dec_max_idx]; 
    1040     } 
    1041 } 
    1042  
    1043 int predict(const model *model_, const feature_node *x) 
    1044 { 
    1045     double *dec_values = Malloc(double, model_->nr_class); 
    1046     int label=predict_values(model_, x, dec_values); 
    1047     free(dec_values); 
    1048     return label; 
    1049 } 
    1050  
    1051 int predict_probability(const struct model *model_, const struct feature_node *x, double* prob_estimates) 
    1052 { 
    1053     if(model_->param.solver_type==L2_LR) 
    1054     { 
    1055         int i; 
    1056         int nr_class=model_->nr_class; 
    1057         int nr_classifier; 
    1058         if(nr_class==2) 
    1059             nr_classifier = 1; 
    1060         else 
    1061             nr_classifier = nr_class; 
    1062  
    1063         int label=predict_values(model_, x, prob_estimates); 
    1064         for(i=0;i<nr_classifier;i++) 
    1065             prob_estimates[i]=1/(1+exp(-prob_estimates[i])); 
    1066  
    1067         if(nr_class==2) // for binary classification 
    1068             prob_estimates[1]=1.-prob_estimates[0]; 
    1069         else 
    1070         { 
    1071             double sum=0; 
    1072             for(i=0; i<nr_class; i++) 
    1073                 sum+=prob_estimates[i]; 
    1074  
    1075             for(i=0; i<nr_class; i++) 
    1076                 prob_estimates[i]=prob_estimates[i]/sum; 
    1077         } 
    1078  
    1079         return label; 
    1080     } 
    1081     else 
    1082         return 0; 
     2270int get_nr_feature(const model *model_) 
     2271{ 
     2272    return model_->nr_feature; 
     2273} 
     2274 
     2275int get_nr_class(const model *model_) 
     2276{ 
     2277    return model_->nr_class; 
     2278} 
     2279 
     2280void get_labels(const model *model_, int* label) 
     2281{ 
     2282    if (model_->label != NULL) 
     2283        for(int i=0;i<model_->nr_class;i++) 
     2284            label[i] = model_->label[i]; 
     2285} 
     2286 
     2287void free_model_content(struct model *model_ptr) 
     2288{ 
     2289    if(model_ptr->w != NULL) 
     2290        free(model_ptr->w); 
     2291    if(model_ptr->label != NULL) 
     2292        free(model_ptr->label); 
     2293} 
     2294 
     2295void free_and_destroy_model(struct model **model_ptr_ptr) 
     2296{ 
     2297    struct model *model_ptr = *model_ptr_ptr; 
     2298    if(model_ptr != NULL) 
     2299    { 
     2300        free_model_content(model_ptr); 
     2301        free(model_ptr); 
     2302    } 
    10832303} 
    10842304 
     
    10952315    if(param->eps <= 0) 
    10962316        return "eps <= 0"; 
    1097  
     2317     
    10982318    if(param->C <= 0) 
    10992319        return "C <= 0"; 
    1100  
    1101     if(param->solver_type != L2_LR 
    1102        && param->solver_type != L2LOSS_SVM_DUAL 
    1103        && param->solver_type != L2LOSS_SVM 
    1104        && param->solver_type != L1LOSS_SVM_DUAL) 
     2320     
     2321    if(param->solver_type != L2R_LR 
     2322       && param->solver_type != L2R_L2LOSS_SVC_DUAL 
     2323       && param->solver_type != L2R_L2LOSS_SVC 
     2324       && param->solver_type != L2R_L1LOSS_SVC_DUAL 
     2325       && param->solver_type != MCSVM_CS 
     2326       && param->solver_type != L1R_L2LOSS_SVC 
     2327       && param->solver_type != L1R_LR 
     2328       && param->solver_type != L2R_LR_DUAL) 
    11052329        return "unknown solver type"; 
    1106  
    1107 //  if(param->solver_type == L1_LR) 
    1108 //      return "sorry! sover_type = 1 (L1_LR) is not supported yet"; 
    1109  
     2330     
    11102331    return NULL; 
    11112332} 
    11122333 
    1113 void cross_validation(const problem *prob, const parameter *param, int nr_fold, int *target) 
    1114 { 
    1115     int i; 
    1116     int *fold_start = Malloc(int,nr_fold+1); 
    1117     int l = prob->l; 
    1118     int *perm = Malloc(int,l); 
    1119  
    1120     for(i=0;i<l;i++) perm[i]=i; 
    1121     for(i=0;i<l;i++) 
    1122     { 
    1123         int j = i+rand()%(l-i); 
    1124         swap(perm[i],perm[j]); 
    1125     } 
    1126     for(i=0;i<=nr_fold;i++) 
    1127         fold_start[i]=i*l/nr_fold; 
    1128  
    1129     for(i=0;i<nr_fold;i++) 
    1130     { 
    1131         int begin = fold_start[i]; 
    1132         int end = fold_start[i+1]; 
    1133         int j,k; 
    1134         struct problem subprob; 
    1135  
    1136         subprob.bias = prob->bias; 
    1137         subprob.n = prob->n; 
    1138         subprob.l = l-(end-begin); 
    1139         subprob.x = Malloc(struct feature_node*,subprob.l); 
    1140         subprob.y = Malloc(int,subprob.l); 
    1141  
    1142         k=0; 
    1143         for(j=0;j<begin;j++) 
    1144         { 
    1145             subprob.x[k] = prob->x[perm[j]]; 
    1146             subprob.y[k] = prob->y[perm[j]]; 
    1147             ++k; 
    1148         } 
    1149         for(j=end;j<l;j++) 
    1150         { 
    1151             subprob.x[k] = prob->x[perm[j]]; 
    1152             subprob.y[k] = prob->y[perm[j]]; 
    1153             ++k; 
    1154         } 
    1155         struct model *submodel = train(&subprob,param); 
    1156         for(j=begin;j<end;j++) 
    1157             target[perm[j]] = predict(submodel,prob->x[perm[j]]); 
    1158         destroy_model(submodel); 
    1159         free(subprob.x); 
    1160         free(subprob.y); 
    1161     } 
    1162     free(fold_start); 
    1163     free(perm); 
    1164 } 
    1165  
    1166 int get_nr_feature(const model *model_) 
    1167 { 
    1168     return model_->nr_feature; 
    1169 } 
    1170  
    1171 int get_nr_class(const model *model_) 
    1172 { 
    1173     return model_->nr_class; 
    1174 } 
    1175  
    1176 void get_labels(const model *model_, int* label) 
    1177 { 
    1178     if (model_->label != NULL) 
    1179         for(int i=0;i<model_->nr_class;i++) 
    1180             label[i] = model_->label[i]; 
    1181 } 
     2334int check_probability_model(const struct model *model_) 
     2335{ 
     2336    return (model_->param.solver_type==L2R_LR || 
     2337            model_->param.solver_type==L2R_LR_DUAL || 
     2338            model_->param.solver_type==L1R_LR); 
     2339} 
     2340 
     2341void set_print_string_function(void (*print_func)(const char*)) 
     2342{ 
     2343    if (print_func == NULL)  
     2344        liblinear_print_string = &print_string_stdout; 
     2345    else 
     2346        liblinear_print_string = print_func; 
     2347} 
     2348 
    11822349 
    11832350/* 
     
    11962363#if _MSC_VER!=0 && _MSC_VER<1300 
    11972364#ifndef min 
    1198 template <class T> inline T min(T x,T y) { return (x<y)?x:y; } 
     2365template <class T> static inline T min(T x,T y) { return (x<y)?x:y; } 
    11992366#endif 
    12002367 
    12012368#ifndef max 
    1202 template <class T> inline T max(T x,T y) { return (x>y)?x:y; } 
     2369template <class T> static inline T max(T x,T y) { return (x>y)?x:y; } 
    12032370#endif 
    12042371#endif 
     
    12072374extern "C" { 
    12082375#endif 
    1209  
    1210 extern double dnrm2_(int *, double *, int *); 
    1211 extern double ddot_(int *, double *, int *, double *, int *); 
    1212 extern int daxpy_(int *, double *, double *, int *, double *, int *); 
    1213 extern int dscal_(int *, double *, double *, int *); 
    1214  
     2376     
     2377    extern double dnrm2_(int *, double *, int *); 
     2378    extern double ddot_(int *, double *, int *, double *, int *); 
     2379    extern int daxpy_(int *, double *, double *, int *, double *, int *); 
     2380    extern int dscal_(int *, double *, double *, int *); 
     2381     
    12152382#ifdef __cplusplus 
    12162383} 
    12172384#endif 
    12182385 
    1219  
    1220 TRON::TRON(const function1 *fun_obj, double eps, int max_iter) 
    1221 { 
    1222     this->fun_obj=const_cast<function1 *>(fun_obj); 
     2386static void default_print(const char *buf) 
     2387{ 
     2388    fputs(buf,stdout); 
     2389    fflush(stdout); 
     2390} 
     2391 
     2392void TRON::info(const char *fmt,...) 
     2393{ 
     2394    char buf[BUFSIZ]; 
     2395    va_list ap; 
     2396    va_start(ap,fmt); 
     2397    vsprintf(buf,fmt,ap); 
     2398    va_end(ap); 
     2399    (*tron_print_string)(buf); 
     2400} 
     2401 
     2402TRON::TRON(const function *fun_obj, double eps, int max_iter) 
     2403{ 
     2404    this->fun_obj=const_cast<function *>(fun_obj); 
    12232405    this->eps=eps; 
    12242406    this->max_iter=max_iter; 
     2407    tron_print_string = default_print; 
    12252408} 
    12262409 
     
    12332416    // Parameters for updating the iterates. 
    12342417    double eta0 = 1e-4, eta1 = 0.25, eta2 = 0.75; 
    1235  
     2418     
    12362419    // Parameters for updating the trust region size delta. 
    12372420    double sigma1 = 0.25, sigma2 = 0.5, sigma3 = 4; 
    1238  
     2421     
    12392422    int n = fun_obj->get_nr_variable(); 
    12402423    int i, cg_iter; 
     
    12462429    double *w_new = new double[n]; 
    12472430    double *g = new double[n]; 
    1248  
     2431     
    12492432    for (i=0; i<n; i++) 
    12502433        w[i] = 0; 
    1251  
    1252         f = fun_obj->fun(w); 
     2434     
     2435    f = fun_obj->fun(w); 
    12532436    fun_obj->grad(w, g); 
    12542437    delta = dnrm2_(&n, g, &inc); 
    12552438    double gnorm1 = delta; 
    12562439    double gnorm = gnorm1; 
    1257  
     2440     
    12582441    if (gnorm <= eps*gnorm1) 
    12592442        search = 0; 
    1260  
     2443     
    12612444    iter = 1; 
    1262  
     2445     
    12632446    while (iter <= max_iter && search) 
    12642447    { 
    12652448        cg_iter = trcg(delta, g, s, r); 
    1266  
     2449         
    12672450        memcpy(w_new, w, sizeof(double)*n); 
    12682451        daxpy_(&n, &one, s, &inc, w_new, &inc); 
    1269  
     2452         
    12702453        gs = ddot_(&n, g, &inc, s, &inc); 
    12712454        prered = -0.5*(gs-ddot_(&n, s, &inc, r, &inc)); 
    1272                 fnew = fun_obj->fun(w_new); 
    1273  
     2455        fnew = fun_obj->fun(w_new); 
     2456         
    12742457        // Compute the actual reduction. 
    1275             actred = f - fnew; 
    1276  
     2458        actred = f - fnew; 
     2459         
    12772460        // On the first iteration, adjust the initial step bound. 
    12782461        snorm = dnrm2_(&n, s, &inc); 
    12792462        if (iter == 1) 
    12802463            delta = min(delta, snorm); 
    1281  
     2464         
    12822465        // Compute prediction alpha*snorm of the step. 
    12832466        if (fnew - f - gs <= 0) 
     
    12852468        else 
    12862469            alpha = max(sigma1, -0.5*(gs/(fnew - f - gs))); 
    1287  
     2470         
    12882471        // Update the trust region bound according to the ratio of actual to predicted reduction. 
    12892472        if (actred < eta0*prered) 
     
    12952478        else 
    12962479            delta = max(delta, min(alpha*snorm, sigma3*delta)); 
    1297  
    1298         //printf("iter %2d act %5.3e pre %5.3e delta %5.3e f %5.3e |g| %5.3e CG %3d\n", iter, actred, prered, delta, f, gnorm, cg_iter); 
    1299  
     2480         
     2481        info("iter %2d act %5.3e pre %5.3e delta %5.3e f %5.3e |g| %5.3e CG %3d\n", iter, actred, prered, delta, f, gnorm, cg_iter); 
     2482         
    13002483        if (actred > eta0*prered) 
    13012484        { 
     
    13032486            memcpy(w, w_new, sizeof(double)*n); 
    13042487            f = fnew; 
    1305                 fun_obj->grad(w, g); 
    1306  
     2488            fun_obj->grad(w, g); 
     2489             
    13072490            gnorm = dnrm2_(&n, g, &inc); 
    13082491            if (gnorm <= eps*gnorm1) 
     
    13112494        if (f < -1.0e+32) 
    13122495        { 
    1313             printf("warning: f < -1.0e+32\n"); 
     2496            info("warning: f < -1.0e+32\n"); 
    13142497            break; 
    13152498        } 
    13162499        if (fabs(actred) <= 0 && prered <= 0) 
    13172500        { 
    1318             printf("warning: actred and prered <= 0\n"); 
     2501            info("warning: actred and prered <= 0\n"); 
    13192502            break; 
    13202503        } 
     
    13222505            fabs(prered) <= 1.0e-12*fabs(f)) 
    13232506        { 
    1324             printf("warning: actred and prered too small\n"); 
     2507            info("warning: actred and prered too small\n"); 
    13252508            break; 
    13262509        } 
    13272510    } 
    1328  
     2511     
    13292512    delete[] g; 
    13302513    delete[] r; 
     
    13412524    double *Hd = new double[n]; 
    13422525    double rTr, rnewTrnew, alpha, beta, cgtol; 
    1343  
     2526     
    13442527    for (i=0; i<n; i++) 
    13452528    { 
     
    13492532    } 
    13502533    cgtol = 0.1*dnrm2_(&n, g, &inc); 
    1351  
     2534     
    13522535    int cg_iter = 0; 
    13532536    rTr = ddot_(&n, r, &inc, r, &inc); 
     
    13582541        cg_iter++; 
    13592542        fun_obj->Hv(d, Hd); 
    1360  
     2543         
    13612544        alpha = rTr/ddot_(&n, d, &inc, Hd, &inc); 
    13622545        daxpy_(&n, &alpha, d, &inc, s, &inc); 
    13632546        if (dnrm2_(&n, s, &inc) > delta) 
    13642547        { 
    1365             //printf("cg reaches trust region boundary\n"); 
     2548            info("cg reaches trust region boundary\n"); 
    13662549            alpha = -alpha; 
    13672550            daxpy_(&n, &alpha, d, &inc, s, &inc); 
    1368  
     2551             
    13692552            double std = ddot_(&n, s, &inc, d, &inc); 
    13702553            double sts = ddot_(&n, s, &inc, s, &inc); 
     
    13892572        rTr = rnewTrnew; 
    13902573    } 
    1391  
     2574     
    13922575    delete[] d; 
    13932576    delete[] Hd; 
    1394  
     2577     
    13952578    return(cg_iter); 
    13962579} 
     
    14052588} 
    14062589 
     2590void TRON::set_print_string(void (*print_string) (const char *buf)) 
     2591{ 
     2592    tron_print_string = print_string; 
     2593} 
     2594 
    14072595 
    14082596/* 
    1409 The folowing load save function1s are used for orange pickling 
     2597The folowing load save functions are used for orange pickling 
    14102598*/ 
    14112599 
     
    14252613 
    14262614    int nr_classifier; 
    1427     if(model_->nr_class==2) 
     2615    if(model_->nr_class==2 && model_->param.solver_type != MCSVM_CS) 
    14282616        nr_classifier=1; 
    14292617    else 
     
    15432731 
    15442732    int nr_classifier; 
    1545     if(nr_class==2) 
     2733    if(nr_class==2 && param.solver_type != MCSVM_CS) 
    15462734        nr_classifier = 1; 
    15472735    else 
     
    16592847} 
    16602848 
     2849static void dont_print_string(const char *s){} 
     2850 
    16612851TLinearLearner::TLinearLearner(){ 
    1662     solver_type = L2_LR; 
     2852    solver_type = L2R_LR; 
    16632853    eps = 0.01f; 
    16642854    C=1; 
     2855    set_print_string_function(&dont_print_string); 
    16652856} 
    16662857 
     
    16972888    domain = examples->domain; 
    16982889    indexMap = _indexMap; 
    1699     computesProbabilities = linmodel->param.solver_type == L2_LR; 
    1700     int nr_classifier = (linmodel->nr_class==2)? 1 : linmodel->nr_class; 
     2890    computesProbabilities = check_probability_model(linmodel); 
     2891    int nr_classifier = (linmodel->nr_class==2 && linmodel->param.solver_type != MCSVM_CS)? 1 : linmodel->nr_class; 
    17012892    weights = mlnew TFloatListList(nr_classifier); 
    17022893    for (int i=0; i<nr_classifier; i++){ 
     
    17092900TLinearClassifier::~TLinearClassifier(){ 
    17102901    if (linmodel) 
    1711         destroy_model(linmodel); 
     2902        free_and_destroy_model(&linmodel); 
    17122903    if (indexMap) 
    17132904        delete indexMap; 
  • source/orange/linear.hpp

    r6796 r7812  
    11/* 
    2 Copyright (c) 2007-2008 The LIBLINEAR Project. 
    3 All rights reserved. 
    4  
    5 Redistribution and use in source and binary forms, with or without 
    6 modification, are permitted provided that the following conditions 
    7 are met: 
    8  
    9 1. Redistributions of source code must retain the above copyright 
    10 notice, this list of conditions and the following disclaimer. 
    11  
    12 2. Redistributions in binary form must reproduce the above copyright 
    13 notice, this list of conditions and the following disclaimer in the 
    14 documentation and/or other materials provided with the distribution. 
    15  
    16 3. Neither name of copyright holders nor the names of its contributors 
    17 may be used to endorse or promote products derived from this software 
    18 without specific prior written permission. 
    19  
    20  
    21 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 
    22 ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 
    23 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR 
    24 A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR 
    25 CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 
    26 EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 
    27 PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 
    28 PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 
    29 LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 
    30 NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 
    31 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 
     2  
     3 Copyright (c) 2007-2010 The LIBLINEAR Project. 
     4 All rights reserved. 
     5  
     6 Redistribution and use in source and binary forms, with or without 
     7 modification, are permitted provided that the following conditions 
     8 are met: 
     9  
     10 1. Redistributions of source code must retain the above copyright 
     11 notice, this list of conditions and the following disclaimer. 
     12  
     13 2. Redistributions in binary form must reproduce the above copyright 
     14 notice, this list of conditions and the following disclaimer in the 
     15 documentation and/or other materials provided with the distribution. 
     16  
     17 3. Neither name of copyright holders nor the names of its contributors 
     18 may be used to endorse or promote products derived from this software 
     19 without specific prior written permission. 
     20  
     21  
     22 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 
     23 ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 
     24 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR 
     25 A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR 
     26 CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 
     27 EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 
     28 PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 
     29 PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 
     30 LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 
     31 NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 
     32 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 
     33  
    3234*/ 
    3335 
     
    4042extern "C" { 
    4143#endif 
    42  
    43 struct feature_node 
    44 { 
    45     int index; 
    46     double value; 
    47 }; 
    48  
    49 struct problem 
    50 { 
    51     int l, n; 
    52     int *y; 
    53     struct feature_node **x; 
    54     double bias;            /* < 0 if no bias term */   
    55 }; 
    56  
    57 enum { L2_LR, L2LOSS_SVM_DUAL, L2LOSS_SVM, L1LOSS_SVM_DUAL }; /* solver_type */ 
    58  
    59 struct parameter 
    60 { 
    61     int solver_type; 
    62  
    63     /* these are for training only */ 
    64     double eps;         /* stopping criteria */ 
    65     double C; 
    66     int nr_weight; 
    67     int *weight_label; 
    68     double* weight; 
    69 }; 
    70  
    71 struct model 
    72 { 
    73     struct parameter param; 
    74     int nr_class;       /* number of classes */ 
    75     int nr_feature; 
    76     double *w; 
    77     int *label;     /* label of each class (label[n]) */ 
    78     double bias; 
    79 }; 
    80  
    81 struct model* train(const struct problem *prob, const struct parameter *param); 
    82 void cross_validation(const struct problem *prob, const struct parameter *param, int nr_fold, int *target); 
    83  
    84 int predict_values(const struct model *model_, const struct feature_node *x, double* dec_values); 
    85 int predict(const struct model *model_, const struct feature_node *x); 
    86 int predict_probability(const struct model *model_, const struct feature_node *x, double* prob_estimates); 
    87  
    88 int save_model(const char *model_file_name, const struct model *model_); 
    89 struct model *load_model(const char *model_file_name); 
    90  
    91 int get_nr_feature(const struct model *model_); 
    92 int get_nr_class(const struct model *model_); 
    93 void get_labels(const struct model *model_, int* label); 
    94  
    95 void destroy_model(struct model *model_); 
    96 void destroy_param(struct parameter *param); 
    97 const char *check_parameter(const struct problem *prob, const struct parameter *param); 
    98  
     44     
     45    struct feature_node 
     46    { 
     47        int index; 
     48        double value; 
     49    }; 
     50     
     51    struct problem 
     52    { 
     53        int l, n; 
     54        int *y; 
     55        struct feature_node **x; 
     56        double bias;            /* < 0 if no bias term */   
     57    }; 
     58     
     59    enum { L2R_LR, L2R_L2LOSS_SVC_DUAL, L2R_L2LOSS_SVC, L2R_L1LOSS_SVC_DUAL, MCSVM_CS, L1R_L2LOSS_SVC, L1R_LR, L2R_LR_DUAL }; /* solver_type */ 
     60     
     61    struct parameter 
     62    { 
     63        int solver_type; 
     64         
     65        /* these are for training only */ 
     66        double eps;         /* stopping criteria */ 
     67        double C; 
     68        int nr_weight; 
     69        int *weight_label; 
     70        double* weight; 
     71    }; 
     72     
     73    struct model 
     74    { 
     75        struct parameter param; 
     76        int nr_class;       /* number of classes */ 
     77        int nr_feature; 
     78        double *w; 
     79        int *label;     /* label of each class */ 
     80        double bias; 
     81    }; 
     82     
     83    struct model* train(const struct problem *prob, const struct parameter *param); 
     84    void cross_validation(const struct problem *prob, const struct parameter *param, int nr_fold, int *target); 
     85     
     86    int predict_values(const struct model *model_, const struct feature_node *x, double* dec_values); 
     87    int predict(const struct model *model_, const struct feature_node *x); 
     88    int predict_probability(const struct model *model_, const struct feature_node *x, double* prob_estimates); 
     89     
     90    int save_model(const char *model_file_name, const struct model *model_); 
     91    struct model *load_model(const char *model_file_name); 
     92     
     93    int get_nr_feature(const struct model *model_); 
     94    int get_nr_class(const struct model *model_); 
     95    void get_labels(const struct model *model_, int* label); 
     96     
     97    void free_model_content(struct model *model_ptr); 
     98    void free_and_destroy_model(struct model **model_ptr_ptr); 
     99    void destroy_param(struct parameter *param); 
     100     
     101    const char *check_parameter(const struct problem *prob, const struct parameter *param); 
     102    int check_probability_model(const struct model *model); 
     103    void set_print_string_function(void (*print_func) (const char*)); 
     104     
    99105#ifdef __cplusplus 
    100106} 
     
    106112#define _TRON_H 
    107113 
    108 class function1 
     114class function 
    109115{ 
    110116public: 
     
    112118    virtual void grad(double *w, double *g) = 0 ; 
    113119    virtual void Hv(double *s, double *Hs) = 0 ; 
    114  
     120     
    115121    virtual int get_nr_variable(void) = 0 ; 
    116     virtual ~function1(void){} 
     122    virtual ~function(void){} 
    117123}; 
    118124 
     
    120126{ 
    121127public: 
    122     TRON(const function1 *fun_obj, double eps = 0.1, int max_iter = 1000); 
     128    TRON(const function *fun_obj, double eps = 0.1, int max_iter = 1000); 
    123129    ~TRON(); 
    124  
     130     
    125131    void tron(double *w); 
    126  
     132    void set_print_string(void (*i_print) (const char *buf)); 
     133     
    127134private: 
    128135    int trcg(double delta, double *g, double *s, double *r); 
    129136    double norm_inf(int n, double *x); 
    130  
     137     
    131138    double eps; 
    132139    int max_iter; 
    133     function1 *fun_obj; 
    134 }; 
    135  
    136 #endif /* _TRON_H */ 
     140    function *fun_obj; 
     141    void info(const char *fmt,...); 
     142    void (*tron_print_string)(const char *buf); 
     143}; 
     144#endif 
     145 
    137146 
    138147#ifndef LINEAR_HPP 
     
    157166    __REGISTER_CLASS 
    158167     
    159     CLASSCONSTANTS(Lossfunction1) enum {L2_LR, L2Loss_SVM_Dual, L2Loss_SVM, L1Loss_SVM_Dual }; 
     168    CLASSCONSTANTS(Lossfunction1_old_) enum {L2_LR, L2Loss_SVM_Dual, L2Loss_SVM, L1Loss_SVM_Dual }; //For backwards compatibility with 1.4 version. 
     169    CLASSCONSTANTS(Lossfunction1) enum {L2R_LR, L2R_L2Loss_SVC_Dual, L2R_L2Loss_SVC, L2R_L1Loss_SVC_Dual, MCSVM_CS, L1R_L2Loss_SVC, L1R_LR, L2R_LR_Dual}; 
     170    CLASSCONSTANTS(LIBLINEAR_VERSION: VERSION=170) 
    160171     
    161172    int solver_type;    //P(&LinearLearner_Lossfunction1) Solver type (loss function1) 
Note: See TracChangeset for help on using the changeset viewer.