SVM_libsvm.cpp

Go to the documentation of this file.
00001 /*
00002  * Copyright (c) 2000-2008 Chih-Chung Chang and Chih-Jen Lin
00003  * All rights reserved.
00004  * 
00005  * Redistribution and use in source and binary forms, with or without
00006  * modification, are permitted provided that the following conditions
00007  * are met:
00008  * 
00009  * 1. Redistributions of source code must retain the above copyright
00010  * notice, this list of conditions and the following disclaimer.
00011  * 
00012  * 2. Redistributions in binary form must reproduce the above copyright
00013  * notice, this list of conditions and the following disclaimer in the
00014  * documentation and/or other materials provided with the distribution.
00015  * 
00016  * 3. Neither name of copyright holders nor the names of its contributors
00017  * may be used to endorse or promote products derived from this software
00018  * without specific prior written permission.
00019  * 
00020  * 
00021  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
00022  * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
00023  * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
00024  * A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR
00025  * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00026  * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00027  * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00028  * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00029  * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00030  * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00031  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00032  *
00033  * Shogun specific adjustments (w) 2006-2008 Soeren Sonnenburg
00034  */
00035 
00036 #include "lib/memory.h"
00037 #include "classifier/svm/SVM_libsvm.h"
00038 #include "kernel/Kernel.h"
00039 #include "lib/io.h"
00040 #include "lib/common.h"
00041 #include "lib/Mathematics.h"
00042 
00043 #include <stdio.h>
00044 #include <stdlib.h>
00045 #include <ctype.h>
00046 #include <float.h>
00047 #include <string.h>
00048 #include <stdarg.h>
00049 typedef KERNELCACHE_ELEM Qfloat;
00050 typedef float64_t schar;
00051 #ifndef min
00052 template <class T> inline T min(T x,T y) { return (x<y)?x:y; }
00053 #endif
00054 #ifndef max
00055 template <class T> inline T max(T x,T y) { return (x>y)?x:y; }
00056 #endif
00057 template <class T> inline void swap(T& x, T& y) { T t=x; x=y; y=t; }
00058 template <class S, class T> inline void clone(T*& dst, S* src, int32_t n)
00059 {
00060     dst = new T[n];
00061     memcpy((void *)dst,(void *)src,sizeof(T)*n);
00062 }
00063 #define INF HUGE_VAL
00064 #define TAU 1e-12
00065 #define Malloc(type,n) (type *)malloc((n)*sizeof(type))
00066 
00067 //
00068 // Kernel Cache
00069 //
00070 // l is the number of total data items
00071 // size is the cache size limit in bytes
00072 //
00073 class Cache
00074 {
00075 public:
00076     Cache(int32_t l, int64_t size);
00077     ~Cache();
00078 
00079     // request data [0,len)
00080     // return some position p where [p,len) need to be filled
00081     // (p >= len if nothing needs to be filled)
00082     int32_t get_data(const int32_t index, Qfloat **data, int32_t len);
00083     void swap_index(int32_t i, int32_t j);  // future_option
00084 
00085 private:
00086     int32_t l;
00087     int64_t size;
00088     struct head_t
00089     {
00090         head_t *prev, *next;    // a cicular list
00091         Qfloat *data;
00092         int32_t len;        // data[0,len) is cached in this entry
00093     };
00094 
00095     head_t *head;
00096     head_t lru_head;
00097     void lru_delete(head_t *h);
00098     void lru_insert(head_t *h);
00099 };
00100 
00101 Cache::Cache(int32_t l_, int64_t size_):l(l_),size(size_)
00102 {
00103     head = (head_t *)calloc(l,sizeof(head_t));  // initialized to 0
00104     size /= sizeof(Qfloat);
00105     size -= l * sizeof(head_t) / sizeof(Qfloat);
00106     size = max(size, (int64_t) 2*l);    // cache must be large enough for two columns
00107     lru_head.next = lru_head.prev = &lru_head;
00108 }
00109 
00110 Cache::~Cache()
00111 {
00112     for(head_t *h = lru_head.next; h != &lru_head; h=h->next)
00113         free(h->data);
00114     free(head);
00115 }
00116 
00117 void Cache::lru_delete(head_t *h)
00118 {
00119     // delete from current location
00120     h->prev->next = h->next;
00121     h->next->prev = h->prev;
00122 }
00123 
00124 void Cache::lru_insert(head_t *h)
00125 {
00126     // insert to last position
00127     h->next = &lru_head;
00128     h->prev = lru_head.prev;
00129     h->prev->next = h;
00130     h->next->prev = h;
00131 }
00132 
00133 int32_t Cache::get_data(const int32_t index, Qfloat **data, int32_t len)
00134 {
00135     head_t *h = &head[index];
00136     if(h->len) lru_delete(h);
00137     int32_t more = len - h->len;
00138 
00139     if(more > 0)
00140     {
00141         // free old space
00142         while(size < more)
00143         {
00144             head_t *old = lru_head.next;
00145             lru_delete(old);
00146             free(old->data);
00147             size += old->len;
00148             old->data = 0;
00149             old->len = 0;
00150         }
00151 
00152         // allocate new space
00153         h->data = (Qfloat *)realloc(h->data,sizeof(Qfloat)*len);
00154         size -= more;
00155         swap(h->len,len);
00156     }
00157 
00158     lru_insert(h);
00159     *data = h->data;
00160     return len;
00161 }
00162 
00163 void Cache::swap_index(int32_t i, int32_t j)
00164 {
00165     if(i==j) return;
00166 
00167     if(head[i].len) lru_delete(&head[i]);
00168     if(head[j].len) lru_delete(&head[j]);
00169     swap(head[i].data,head[j].data);
00170     swap(head[i].len,head[j].len);
00171     if(head[i].len) lru_insert(&head[i]);
00172     if(head[j].len) lru_insert(&head[j]);
00173 
00174     if(i>j) swap(i,j);
00175     for(head_t *h = lru_head.next; h!=&lru_head; h=h->next)
00176     {
00177         if(h->len > i)
00178         {
00179             if(h->len > j)
00180                 swap(h->data[i],h->data[j]);
00181             else
00182             {
00183                 // give up
00184                 lru_delete(h);
00185                 free(h->data);
00186                 size += h->len;
00187                 h->data = 0;
00188                 h->len = 0;
00189             }
00190         }
00191     }
00192 }
00193 
00194 //
00195 // Kernel evaluation
00196 //
00197 // the static method k_function is for doing single kernel evaluation
00198 // the constructor of Kernel prepares to calculate the l*l kernel matrix
00199 // the member function get_Q is for getting one column from the Q Matrix
00200 //
00201 class QMatrix {
00202 public:
00203     virtual Qfloat *get_Q(int32_t column, int32_t len) const = 0;
00204     virtual Qfloat *get_QD() const = 0;
00205     virtual void swap_index(int32_t i, int32_t j) const = 0;
00206     virtual ~QMatrix() {}
00207 };
00208 
00209 class Kernel: public QMatrix {
00210 public:
00211     Kernel(int32_t l, svm_node * const * x, const svm_parameter& param);
00212     virtual ~Kernel();
00213 
00214     virtual Qfloat *get_Q(int32_t column, int32_t len) const = 0;
00215     virtual Qfloat *get_QD() const = 0;
00216     virtual void swap_index(int32_t i, int32_t j) const // no so const...
00217     {
00218         swap(x[i],x[j]);
00219         if(x_square) swap(x_square[i],x_square[j]);
00220     }
00221 
00222     inline float64_t kernel_function(int32_t i, int32_t j) const
00223     {
00224         return kernel->kernel(x[i]->index,x[j]->index);
00225     }
00226 
00227 private:
00228     CKernel* kernel;
00229     const svm_node **x;
00230     float64_t *x_square;
00231 
00232     // svm_parameter
00233     const int32_t kernel_type;
00234     const int32_t degree;
00235     const float64_t gamma;
00236     const float64_t coef0;
00237 };
00238 
00239 Kernel::Kernel(int32_t l, svm_node * const * x_, const svm_parameter& param)
00240 : kernel_type(param.kernel_type), degree(param.degree),
00241  gamma(param.gamma), coef0(param.coef0)
00242 {
00243     clone(x,x_,l);
00244     x_square = 0;
00245     kernel=param.kernel;
00246 }
00247 
00248 Kernel::~Kernel()
00249 {
00250     delete[] x;
00251     delete[] x_square;
00252 }
00253 
00254 // Generalized SMO+SVMlight algorithm
00255 // Solves:
00256 //
00257 //  min 0.5(\alpha^T Q \alpha) + b^T \alpha
00258 //
00259 //      y^T \alpha = \delta
00260 //      y_i = +1 or -1
00261 //      0 <= alpha_i <= Cp for y_i = 1
00262 //      0 <= alpha_i <= Cn for y_i = -1
00263 //
00264 // Given:
00265 //
00266 //  Q, p, y, Cp, Cn, and an initial feasible point \alpha
00267 //  l is the size of vectors and matrices
00268 //  eps is the stopping tolerance
00269 //
00270 // solution will be put in \alpha, objective value will be put in obj
00271 //
00272 class Solver {
00273 public:
00274     Solver() {};
00275     virtual ~Solver() {};
00276 
00277     struct SolutionInfo {
00278         float64_t obj;
00279         float64_t rho;
00280         float64_t upper_bound_p;
00281         float64_t upper_bound_n;
00282         float64_t r;    // for Solver_NU
00283     };
00284 
00285     void Solve(
00286         int32_t l, const QMatrix& Q, const float64_t *p_, const schar *y_,
00287         float64_t *alpha_, float64_t Cp, float64_t Cn, float64_t eps,
00288         SolutionInfo* si, int32_t shrinking);
00289 
00290 protected:
00291     int32_t active_size;
00292     schar *y;
00293     float64_t *G;       // gradient of objective function
00294     enum { LOWER_BOUND, UPPER_BOUND, FREE };
00295     char *alpha_status; // LOWER_BOUND, UPPER_BOUND, FREE
00296     float64_t *alpha;
00297     const QMatrix *Q;
00298     const Qfloat *QD;
00299     float64_t eps;
00300     float64_t Cp,Cn;
00301     float64_t *p;
00302     int32_t *active_set;
00303     float64_t *G_bar;       // gradient, if we treat free variables as 0
00304     int32_t l;
00305     bool unshrinked;    // XXX
00306 
00307     float64_t get_C(int32_t i)
00308     {
00309         return (y[i] > 0)? Cp : Cn;
00310     }
00311     void update_alpha_status(int32_t i)
00312     {
00313         if(alpha[i] >= get_C(i))
00314             alpha_status[i] = UPPER_BOUND;
00315         else if(alpha[i] <= 0)
00316             alpha_status[i] = LOWER_BOUND;
00317         else alpha_status[i] = FREE;
00318     }
00319     bool is_upper_bound(int32_t i) { return alpha_status[i] == UPPER_BOUND; }
00320     bool is_lower_bound(int32_t i) { return alpha_status[i] == LOWER_BOUND; }
00321     bool is_free(int32_t i) { return alpha_status[i] == FREE; }
00322     void swap_index(int32_t i, int32_t j);
00323     void reconstruct_gradient();
00324     virtual int32_t select_working_set(int32_t &i, int32_t &j, float64_t &gap);
00325     virtual float64_t calculate_rho();
00326     virtual void do_shrinking();
00327 private:
00328     bool be_shrunken(int32_t i, float64_t Gmax1, float64_t Gmax2);
00329 };
00330 
00331 void Solver::swap_index(int32_t i, int32_t j)
00332 {
00333     Q->swap_index(i,j);
00334     swap(y[i],y[j]);
00335     swap(G[i],G[j]);
00336     swap(alpha_status[i],alpha_status[j]);
00337     swap(alpha[i],alpha[j]);
00338     swap(p[i],p[j]);
00339     swap(active_set[i],active_set[j]);
00340     swap(G_bar[i],G_bar[j]);
00341 }
00342 
00343 void Solver::reconstruct_gradient()
00344 {
00345     // reconstruct inactive elements of G from G_bar and free variables
00346 
00347     if(active_size == l) return;
00348 
00349     int32_t i;
00350     for(i=active_size;i<l;i++)
00351         G[i] = G_bar[i] + p[i];
00352     
00353     for(i=0;i<active_size;i++)
00354         if(is_free(i))
00355         {
00356             const Qfloat *Q_i = Q->get_Q(i,l);
00357             float64_t alpha_i = alpha[i];
00358             for(int32_t j=active_size;j<l;j++)
00359                 G[j] += alpha_i * Q_i[j];
00360         }
00361 }
00362 
00363 void Solver::Solve(
00364     int32_t p_l, const QMatrix& p_Q, const float64_t *p_p,
00365     const schar *p_y, float64_t *p_alpha, float64_t p_Cp, float64_t p_Cn,
00366     float64_t p_eps, SolutionInfo* p_si, int32_t shrinking)
00367 {
00368     this->l = p_l;
00369     this->Q = &p_Q;
00370     QD=Q->get_QD();
00371     clone(p, p_p,l);
00372     clone(y, p_y,l);
00373     clone(alpha,p_alpha,l);
00374     this->Cp = p_Cp;
00375     this->Cn = p_Cn;
00376     this->eps = p_eps;
00377     unshrinked = false;
00378 
00379     // initialize alpha_status
00380     {
00381         alpha_status = new char[l];
00382         for(int32_t i=0;i<l;i++)
00383             update_alpha_status(i);
00384     }
00385 
00386     // initialize active set (for shrinking)
00387     {
00388         active_set = new int32_t[l];
00389         for(int32_t i=0;i<l;i++)
00390             active_set[i] = i;
00391         active_size = l;
00392     }
00393 
00394     // initialize gradient
00395     {
00396         G = new float64_t[l];
00397         G_bar = new float64_t[l];
00398         int32_t i;
00399         for(i=0;i<l;i++)
00400         {
00401             G[i] = p_p[i];
00402             G_bar[i] = 0;
00403         }
00404         for(i=0;i<l;i++)
00405             if(!is_lower_bound(i))
00406             {
00407                 const Qfloat *Q_i = Q->get_Q(i,l);
00408                 float64_t alpha_i = alpha[i];
00409                 int32_t j;
00410                 for(j=0;j<l;j++)
00411                     G[j] += alpha_i*Q_i[j];
00412                 if(is_upper_bound(i))
00413                     for(j=0;j<l;j++)
00414                         G_bar[j] += get_C(i) * Q_i[j];
00415             }
00416     }
00417 
00418     // optimization step
00419 
00420     int32_t iter = 0;
00421     int32_t counter = min(l,1000)+1;
00422 
00423     while(1)
00424     {
00425         // show progress and do shrinking
00426 
00427         if(--counter == 0)
00428         {
00429             counter = min(l,1000);
00430             if(shrinking) do_shrinking();
00431             //SG_SINFO(".");
00432         }
00433 
00434         int32_t i,j;
00435         float64_t gap;
00436         if(select_working_set(i,j, gap)!=0)
00437         {
00438             // reconstruct the whole gradient
00439             reconstruct_gradient();
00440             // reset active set size and check
00441             active_size = l;
00442             //SG_SINFO("*");
00443             if(select_working_set(i,j, gap)!=0)
00444                 break;
00445             else
00446                 counter = 1;    // do shrinking next iteration
00447         }
00448 
00449         SG_SABS_PROGRESS(gap, -CMath::log10(gap), -CMath::log10(1), -CMath::log10(eps), 6);
00450         
00451         ++iter;
00452 
00453         // update alpha[i] and alpha[j], handle bounds carefully
00454         
00455         const Qfloat *Q_i = Q->get_Q(i,active_size);
00456         const Qfloat *Q_j = Q->get_Q(j,active_size);
00457 
00458         float64_t C_i = get_C(i);
00459         float64_t C_j = get_C(j);
00460 
00461         float64_t old_alpha_i = alpha[i];
00462         float64_t old_alpha_j = alpha[j];
00463 
00464         if(y[i]!=y[j])
00465         {
00466             float64_t quad_coef = Q_i[i]+Q_j[j]+2*Q_i[j];
00467             if (quad_coef <= 0)
00468                 quad_coef = TAU;
00469             float64_t delta = (-G[i]-G[j])/quad_coef;
00470             float64_t diff = alpha[i] - alpha[j];
00471             alpha[i] += delta;
00472             alpha[j] += delta;
00473             
00474             if(diff > 0)
00475             {
00476                 if(alpha[j] < 0)
00477                 {
00478                     alpha[j] = 0;
00479                     alpha[i] = diff;
00480                 }
00481             }
00482             else
00483             {
00484                 if(alpha[i] < 0)
00485                 {
00486                     alpha[i] = 0;
00487                     alpha[j] = -diff;
00488                 }
00489             }
00490             if(diff > C_i - C_j)
00491             {
00492                 if(alpha[i] > C_i)
00493                 {
00494                     alpha[i] = C_i;
00495                     alpha[j] = C_i - diff;
00496                 }
00497             }
00498             else
00499             {
00500                 if(alpha[j] > C_j)
00501                 {
00502                     alpha[j] = C_j;
00503                     alpha[i] = C_j + diff;
00504                 }
00505             }
00506         }
00507         else
00508         {
00509             float64_t quad_coef = Q_i[i]+Q_j[j]-2*Q_i[j];
00510             if (quad_coef <= 0)
00511                 quad_coef = TAU;
00512             float64_t delta = (G[i]-G[j])/quad_coef;
00513             float64_t sum = alpha[i] + alpha[j];
00514             alpha[i] -= delta;
00515             alpha[j] += delta;
00516 
00517             if(sum > C_i)
00518             {
00519                 if(alpha[i] > C_i)
00520                 {
00521                     alpha[i] = C_i;
00522                     alpha[j] = sum - C_i;
00523                 }
00524             }
00525             else
00526             {
00527                 if(alpha[j] < 0)
00528                 {
00529                     alpha[j] = 0;
00530                     alpha[i] = sum;
00531                 }
00532             }
00533             if(sum > C_j)
00534             {
00535                 if(alpha[j] > C_j)
00536                 {
00537                     alpha[j] = C_j;
00538                     alpha[i] = sum - C_j;
00539                 }
00540             }
00541             else
00542             {
00543                 if(alpha[i] < 0)
00544                 {
00545                     alpha[i] = 0;
00546                     alpha[j] = sum;
00547                 }
00548             }
00549         }
00550 
00551         // update G
00552 
00553         float64_t delta_alpha_i = alpha[i] - old_alpha_i;
00554         float64_t delta_alpha_j = alpha[j] - old_alpha_j;
00555         
00556         for(int32_t k=0;k<active_size;k++)
00557         {
00558             G[k] += Q_i[k]*delta_alpha_i + Q_j[k]*delta_alpha_j;
00559         }
00560 
00561         // update alpha_status and G_bar
00562 
00563         {
00564             bool ui = is_upper_bound(i);
00565             bool uj = is_upper_bound(j);
00566             update_alpha_status(i);
00567             update_alpha_status(j);
00568             int32_t k;
00569             if(ui != is_upper_bound(i))
00570             {
00571                 Q_i = Q->get_Q(i,l);
00572                 if(ui)
00573                     for(k=0;k<l;k++)
00574                         G_bar[k] -= C_i * Q_i[k];
00575                 else
00576                     for(k=0;k<l;k++)
00577                         G_bar[k] += C_i * Q_i[k];
00578             }
00579 
00580             if(uj != is_upper_bound(j))
00581             {
00582                 Q_j = Q->get_Q(j,l);
00583                 if(uj)
00584                     for(k=0;k<l;k++)
00585                         G_bar[k] -= C_j * Q_j[k];
00586                 else
00587                     for(k=0;k<l;k++)
00588                         G_bar[k] += C_j * Q_j[k];
00589             }
00590         }
00591     }
00592 
00593     // calculate rho
00594 
00595     p_si->rho = calculate_rho();
00596 
00597     // calculate objective value
00598     {
00599         float64_t v = 0;
00600         int32_t i;
00601         for(i=0;i<l;i++)
00602             v += alpha[i] * (G[i] + p[i]);
00603 
00604         p_si->obj = v/2;
00605     }
00606 
00607     // put back the solution
00608     {
00609         for(int32_t i=0;i<l;i++)
00610             p_alpha[active_set[i]] = alpha[i];
00611     }
00612 
00613     p_si->upper_bound_p = Cp;
00614     p_si->upper_bound_n = Cn;
00615 
00616     SG_SINFO("\noptimization finished, #iter = %d\n",iter);
00617 
00618     delete[] p;
00619     delete[] y;
00620     delete[] alpha;
00621     delete[] alpha_status;
00622     delete[] active_set;
00623     delete[] G;
00624     delete[] G_bar;
00625 }
00626 
00627 // return 1 if already optimal, return 0 otherwise
00628 int32_t Solver::select_working_set(
00629     int32_t &out_i, int32_t &out_j, float64_t &gap)
00630 {
00631     // return i,j such that
00632     // i: maximizes -y_i * grad(f)_i, i in I_up(\alpha)
00633     // j: minimizes the decrease of obj value
00634     //    (if quadratic coefficient <= 0, replace it with tau)
00635     //    -y_j*grad(f)_j < -y_i*grad(f)_i, j in I_low(\alpha)
00636     
00637     float64_t Gmax = -INF;
00638     float64_t Gmax2 = -INF;
00639     int32_t Gmax_idx = -1;
00640     int32_t Gmin_idx = -1;
00641     float64_t obj_diff_min = INF;
00642 
00643     for(int32_t t=0;t<active_size;t++)
00644         if(y[t]==+1)
00645         {
00646             if(!is_upper_bound(t))
00647                 if(-G[t] >= Gmax)
00648                 {
00649                     Gmax = -G[t];
00650                     Gmax_idx = t;
00651                 }
00652         }
00653         else
00654         {
00655             if(!is_lower_bound(t))
00656                 if(G[t] >= Gmax)
00657                 {
00658                     Gmax = G[t];
00659                     Gmax_idx = t;
00660                 }
00661         }
00662 
00663     int32_t i = Gmax_idx;
00664     const Qfloat *Q_i = NULL;
00665     if(i != -1) // NULL Q_i not accessed: Gmax=-INF if i=-1
00666         Q_i = Q->get_Q(i,active_size);
00667 
00668     for(int32_t j=0;j<active_size;j++)
00669     {
00670         if(y[j]==+1)
00671         {
00672             if (!is_lower_bound(j))
00673             {
00674                 float64_t grad_diff=Gmax+G[j];
00675                 if (G[j] >= Gmax2)
00676                     Gmax2 = G[j];
00677                 if (grad_diff > 0)
00678                 {
00679                     float64_t obj_diff; 
00680                     float64_t quad_coef=Q_i[i]+QD[j]-2*y[i]*Q_i[j];
00681                     if (quad_coef > 0)
00682                         obj_diff = -(grad_diff*grad_diff)/quad_coef;
00683                     else
00684                         obj_diff = -(grad_diff*grad_diff)/TAU;
00685 
00686                     if (obj_diff <= obj_diff_min)
00687                     {
00688                         Gmin_idx=j;
00689                         obj_diff_min = obj_diff;
00690                     }
00691                 }
00692             }
00693         }
00694         else
00695         {
00696             if (!is_upper_bound(j))
00697             {
00698                 float64_t grad_diff= Gmax-G[j];
00699                 if (-G[j] >= Gmax2)
00700                     Gmax2 = -G[j];
00701                 if (grad_diff > 0)
00702                 {
00703                     float64_t obj_diff; 
00704                     float64_t quad_coef=Q_i[i]+QD[j]+2*y[i]*Q_i[j];
00705                     if (quad_coef > 0)
00706                         obj_diff = -(grad_diff*grad_diff)/quad_coef;
00707                     else
00708                         obj_diff = -(grad_diff*grad_diff)/TAU;
00709 
00710                     if (obj_diff <= obj_diff_min)
00711                     {
00712                         Gmin_idx=j;
00713                         obj_diff_min = obj_diff;
00714                     }
00715                 }
00716             }
00717         }
00718     }
00719 
00720     gap=Gmax+Gmax2;
00721     if(gap < eps)
00722         return 1;
00723 
00724     out_i = Gmax_idx;
00725     out_j = Gmin_idx;
00726     return 0;
00727 }
00728 
00729 bool Solver::be_shrunken(int32_t i, float64_t Gmax1, float64_t Gmax2)
00730 {
00731     if(is_upper_bound(i))
00732     {
00733         if(y[i]==+1)
00734             return(-G[i] > Gmax1);
00735         else
00736             return(-G[i] > Gmax2);
00737     }
00738     else if(is_lower_bound(i))
00739     {
00740         if(y[i]==+1)
00741             return(G[i] > Gmax2);
00742         else    
00743             return(G[i] > Gmax1);
00744     }
00745     else
00746         return(false);
00747 }
00748 
00749 void Solver::do_shrinking()
00750 {
00751     int32_t i;
00752     float64_t Gmax1 = -INF;     // max { -y_i * grad(f)_i | i in I_up(\alpha) }
00753     float64_t Gmax2 = -INF;     // max { y_i * grad(f)_i | i in I_low(\alpha) }
00754 
00755     // find maximal violating pair first
00756     for(i=0;i<active_size;i++)
00757     {
00758         if(y[i]==+1)    
00759         {
00760             if(!is_upper_bound(i))  
00761             {
00762                 if(-G[i] >= Gmax1)
00763                     Gmax1 = -G[i];
00764             }
00765             if(!is_lower_bound(i))  
00766             {
00767                 if(G[i] >= Gmax2)
00768                     Gmax2 = G[i];
00769             }
00770         }
00771         else    
00772         {
00773             if(!is_upper_bound(i))  
00774             {
00775                 if(-G[i] >= Gmax2)
00776                     Gmax2 = -G[i];
00777             }
00778             if(!is_lower_bound(i))  
00779             {
00780                 if(G[i] >= Gmax1)
00781                     Gmax1 = G[i];
00782             }
00783         }
00784     }
00785 
00786     // shrink
00787 
00788     for(i=0;i<active_size;i++)
00789         if (be_shrunken(i, Gmax1, Gmax2))
00790         {
00791             active_size--;
00792             while (active_size > i)
00793             {
00794                 if (!be_shrunken(active_size, Gmax1, Gmax2))
00795                 {
00796                     swap_index(i,active_size);
00797                     break;
00798                 }
00799                 active_size--;
00800             }
00801         }
00802 
00803     // unshrink, check all variables again before final iterations
00804 
00805     if(unshrinked || Gmax1 + Gmax2 > eps*10) return;
00806     
00807     unshrinked = true;
00808     reconstruct_gradient();
00809 
00810     for(i=l-1;i>=active_size;i--)
00811         if (!be_shrunken(i, Gmax1, Gmax2))
00812         {
00813             while (active_size < i)
00814             {
00815                 if (be_shrunken(active_size, Gmax1, Gmax2))
00816                 {
00817                     swap_index(i,active_size);
00818                     break;
00819                 }
00820                 active_size++;
00821             }
00822             active_size++;
00823         }
00824 }
00825 
00826 float64_t Solver::calculate_rho()
00827 {
00828     float64_t r;
00829     int32_t nr_free = 0;
00830     float64_t ub = INF, lb = -INF, sum_free = 0;
00831     for(int32_t i=0;i<active_size;i++)
00832     {
00833         float64_t yG = y[i]*G[i];
00834 
00835         if(is_upper_bound(i))
00836         {
00837             if(y[i]==-1)
00838                 ub = min(ub,yG);
00839             else
00840                 lb = max(lb,yG);
00841         }
00842         else if(is_lower_bound(i))
00843         {
00844             if(y[i]==+1)
00845                 ub = min(ub,yG);
00846             else
00847                 lb = max(lb,yG);
00848         }
00849         else
00850         {
00851             ++nr_free;
00852             sum_free += yG;
00853         }
00854     }
00855 
00856     if(nr_free>0)
00857         r = sum_free/nr_free;
00858     else
00859         r = (ub+lb)/2;
00860 
00861     return r;
00862 }
00863 
00864 //
00865 // Solver for nu-svm classification and regression
00866 //
00867 // additional constraint: e^T \alpha = constant
00868 //
00869 class Solver_NU : public Solver
00870 {
00871 public:
00872     Solver_NU() {}
00873     void Solve(
00874         int32_t p_l, const QMatrix& p_Q, const float64_t *p_p,
00875         const schar *p_y, float64_t* p_alpha, float64_t p_Cp, float64_t p_Cn,
00876         float64_t p_eps, SolutionInfo* p_si, int32_t shrinking)
00877     {
00878         this->si = p_si;
00879         Solver::Solve(p_l,p_Q,p,p_y,p_alpha,p_Cp,p_Cn,p_eps,p_si,shrinking);
00880     }
00881 private:
00882     SolutionInfo *si;
00883     int32_t select_working_set(int32_t &i, int32_t &j, float64_t &gap);
00884     float64_t calculate_rho();
00885     bool be_shrunken(
00886         int32_t i, float64_t Gmax1, float64_t Gmax2, float64_t Gmax3,
00887         float64_t Gmax4);
00888     void do_shrinking();
00889 };
00890 
00891 // return 1 if already optimal, return 0 otherwise
00892 int32_t Solver_NU::select_working_set(
00893     int32_t &out_i, int32_t &out_j, float64_t &gap)
00894 {
00895     // return i,j such that y_i = y_j and
00896     // i: maximizes -y_i * grad(f)_i, i in I_up(\alpha)
00897     // j: minimizes the decrease of obj value
00898     //    (if quadratic coefficient <= 0, replace it with tau)
00899     //    -y_j*grad(f)_j < -y_i*grad(f)_i, j in I_low(\alpha)
00900 
00901     float64_t Gmaxp = -INF;
00902     float64_t Gmaxp2 = -INF;
00903     int32_t Gmaxp_idx = -1;
00904 
00905     float64_t Gmaxn = -INF;
00906     float64_t Gmaxn2 = -INF;
00907     int32_t Gmaxn_idx = -1;
00908 
00909     int32_t Gmin_idx = -1;
00910     float64_t obj_diff_min = INF;
00911 
00912     for(int32_t t=0;t<active_size;t++)
00913         if(y[t]==+1)
00914         {
00915             if(!is_upper_bound(t))
00916                 if(-G[t] >= Gmaxp)
00917                 {
00918                     Gmaxp = -G[t];
00919                     Gmaxp_idx = t;
00920                 }
00921         }
00922         else
00923         {
00924             if(!is_lower_bound(t))
00925                 if(G[t] >= Gmaxn)
00926                 {
00927                     Gmaxn = G[t];
00928                     Gmaxn_idx = t;
00929                 }
00930         }
00931 
00932     int32_t ip = Gmaxp_idx;
00933     int32_t in = Gmaxn_idx;
00934     const Qfloat *Q_ip = NULL;
00935     const Qfloat *Q_in = NULL;
00936     if(ip != -1) // NULL Q_ip not accessed: Gmaxp=-INF if ip=-1
00937         Q_ip = Q->get_Q(ip,active_size);
00938     if(in != -1)
00939         Q_in = Q->get_Q(in,active_size);
00940 
00941     for(int32_t j=0;j<active_size;j++)
00942     {
00943         if(y[j]==+1)
00944         {
00945             if (!is_lower_bound(j)) 
00946             {
00947                 float64_t grad_diff=Gmaxp+G[j];
00948                 if (G[j] >= Gmaxp2)
00949                     Gmaxp2 = G[j];
00950                 if (grad_diff > 0)
00951                 {
00952                     float64_t obj_diff; 
00953                     float64_t quad_coef = Q_ip[ip]+QD[j]-2*Q_ip[j];
00954                     if (quad_coef > 0)
00955                         obj_diff = -(grad_diff*grad_diff)/quad_coef;
00956                     else
00957                         obj_diff = -(grad_diff*grad_diff)/TAU;
00958 
00959                     if (obj_diff <= obj_diff_min)
00960                     {
00961                         Gmin_idx=j;
00962                         obj_diff_min = obj_diff;
00963                     }
00964                 }
00965             }
00966         }
00967         else
00968         {
00969             if (!is_upper_bound(j))
00970             {
00971                 float64_t grad_diff=Gmaxn-G[j];
00972                 if (-G[j] >= Gmaxn2)
00973                     Gmaxn2 = -G[j];
00974                 if (grad_diff > 0)
00975                 {
00976                     float64_t obj_diff; 
00977                     float64_t quad_coef = Q_in[in]+QD[j]-2*Q_in[j];
00978                     if (quad_coef > 0)
00979                         obj_diff = -(grad_diff*grad_diff)/quad_coef;
00980                     else
00981                         obj_diff = -(grad_diff*grad_diff)/TAU;
00982 
00983                     if (obj_diff <= obj_diff_min)
00984                     {
00985                         Gmin_idx=j;
00986                         obj_diff_min = obj_diff;
00987                     }
00988                 }
00989             }
00990         }
00991     }
00992 
00993     gap=max(Gmaxp+Gmaxp2,Gmaxn+Gmaxn2);
00994     if(gap < eps)
00995         return 1;
00996 
00997     if (y[Gmin_idx] == +1)
00998         out_i = Gmaxp_idx;
00999     else
01000         out_i = Gmaxn_idx;
01001     out_j = Gmin_idx;
01002 
01003     return 0;
01004 }
01005 
01006 bool Solver_NU::be_shrunken(
01007     int32_t i, float64_t Gmax1, float64_t Gmax2, float64_t Gmax3,
01008     float64_t Gmax4)
01009 {
01010     if(is_upper_bound(i))
01011     {
01012         if(y[i]==+1)
01013             return(-G[i] > Gmax1);
01014         else    
01015             return(-G[i] > Gmax4);
01016     }
01017     else if(is_lower_bound(i))
01018     {
01019         if(y[i]==+1)
01020             return(G[i] > Gmax2);
01021         else    
01022             return(G[i] > Gmax3);
01023     }
01024     else
01025         return(false);
01026 }
01027 
01028 void Solver_NU::do_shrinking()
01029 {
01030     float64_t Gmax1 = -INF; // max { -y_i * grad(f)_i | y_i = +1, i in I_up(\alpha) }
01031     float64_t Gmax2 = -INF; // max { y_i * grad(f)_i | y_i = +1, i in I_low(\alpha) }
01032     float64_t Gmax3 = -INF; // max { -y_i * grad(f)_i | y_i = -1, i in I_up(\alpha) }
01033     float64_t Gmax4 = -INF; // max { y_i * grad(f)_i | y_i = -1, i in I_low(\alpha) }
01034 
01035     // find maximal violating pair first
01036     int32_t i;
01037     for(i=0;i<active_size;i++)
01038     {
01039         if(!is_upper_bound(i))
01040         {
01041             if(y[i]==+1)
01042             {
01043                 if(-G[i] > Gmax1) Gmax1 = -G[i];
01044             }
01045             else    if(-G[i] > Gmax4) Gmax4 = -G[i];
01046         }
01047         if(!is_lower_bound(i))
01048         {
01049             if(y[i]==+1)
01050             {   
01051                 if(G[i] > Gmax2) Gmax2 = G[i];
01052             }
01053             else    if(G[i] > Gmax3) Gmax3 = G[i];
01054         }
01055     }
01056 
01057     // shrinking
01058 
01059     for(i=0;i<active_size;i++)
01060         if (be_shrunken(i, Gmax1, Gmax2, Gmax3, Gmax4))
01061         {
01062             active_size--;
01063             while (active_size > i)
01064             {
01065                 if (!be_shrunken(active_size, Gmax1, Gmax2, Gmax3, Gmax4))
01066                 {
01067                     swap_index(i,active_size);
01068                     break;
01069                 }
01070                 active_size--;
01071             }
01072         }
01073 
01074     // unshrink, check all variables again before final iterations
01075 
01076     if(unshrinked || max(Gmax1+Gmax2,Gmax3+Gmax4) > eps*10) return;
01077     
01078     unshrinked = true;
01079     reconstruct_gradient();
01080 
01081     for(i=l-1;i>=active_size;i--)
01082         if (!be_shrunken(i, Gmax1, Gmax2, Gmax3, Gmax4))
01083         {
01084             while (active_size < i)
01085             {
01086                 if (be_shrunken(active_size, Gmax1, Gmax2, Gmax3, Gmax4))
01087                 {
01088                     swap_index(i,active_size);
01089                     break;
01090                 }
01091                 active_size++;
01092             }
01093             active_size++;
01094         }
01095 }
01096 
01097 float64_t Solver_NU::calculate_rho()
01098 {
01099     int32_t nr_free1 = 0,nr_free2 = 0;
01100     float64_t ub1 = INF, ub2 = INF;
01101     float64_t lb1 = -INF, lb2 = -INF;
01102     float64_t sum_free1 = 0, sum_free2 = 0;
01103 
01104     for(int32_t i=0; i<active_size; i++)
01105     {
01106         if(y[i]==+1)
01107         {
01108             if(is_upper_bound(i))
01109                 lb1 = max(lb1,G[i]);
01110             else if(is_lower_bound(i))
01111                 ub1 = min(ub1,G[i]);
01112             else
01113             {
01114                 ++nr_free1;
01115                 sum_free1 += G[i];
01116             }
01117         }
01118         else
01119         {
01120             if(is_upper_bound(i))
01121                 lb2 = max(lb2,G[i]);
01122             else if(is_lower_bound(i))
01123                 ub2 = min(ub2,G[i]);
01124             else
01125             {
01126                 ++nr_free2;
01127                 sum_free2 += G[i];
01128             }
01129         }
01130     }
01131 
01132     float64_t r1,r2;
01133     if(nr_free1 > 0)
01134         r1 = sum_free1/nr_free1;
01135     else
01136         r1 = (ub1+lb1)/2;
01137     
01138     if(nr_free2 > 0)
01139         r2 = sum_free2/nr_free2;
01140     else
01141         r2 = (ub2+lb2)/2;
01142     
01143     si->r = (r1+r2)/2;
01144     return (r1-r2)/2;
01145 }
01146 
01147 //
01148 // Q matrices for various formulations
01149 //
01150 class SVC_Q: public Kernel
01151 { 
01152 public:
01153     SVC_Q(const svm_problem& prob, const svm_parameter& param, const schar *y_)
01154     :Kernel(prob.l, prob.x, param)
01155     {
01156         clone(y,y_,prob.l);
01157         cache = new Cache(prob.l,(int64_t)(param.cache_size*(1l<<20)));
01158         QD = new Qfloat[prob.l];
01159         for(int32_t i=0;i<prob.l;i++)
01160             QD[i]= (Qfloat)kernel_function(i,i);
01161     }
01162     
01163     Qfloat *get_Q(int32_t i, int32_t len) const
01164     {
01165         Qfloat *data;
01166         int32_t start;
01167         if((start = cache->get_data(i,&data,len)) < len)
01168         {
01169             for(int32_t j=start;j<len;j++)
01170                 data[j] = (Qfloat) y[i]*y[j]*kernel_function(i,j);
01171         }
01172         return data;
01173     }
01174 
01175     Qfloat *get_QD() const
01176     {
01177         return QD;
01178     }
01179 
01180     void swap_index(int32_t i, int32_t j) const
01181     {
01182         cache->swap_index(i,j);
01183         Kernel::swap_index(i,j);
01184         swap(y[i],y[j]);
01185         swap(QD[i],QD[j]);
01186     }
01187 
01188     ~SVC_Q()
01189     {
01190         delete[] y;
01191         delete cache;
01192         delete[] QD;
01193     }
01194 private:
01195     schar *y;
01196     Cache *cache;
01197     Qfloat *QD;
01198 };
01199 
01200 class ONE_CLASS_Q: public Kernel
01201 {
01202 public:
01203     ONE_CLASS_Q(const svm_problem& prob, const svm_parameter& param)
01204     :Kernel(prob.l, prob.x, param)
01205     {
01206         cache = new Cache(prob.l,(int64_t)(param.cache_size*(1l<<20)));
01207         QD = new Qfloat[prob.l];
01208         for(int32_t i=0;i<prob.l;i++)
01209             QD[i]= (Qfloat)kernel_function(i,i);
01210     }
01211     
01212     Qfloat *get_Q(int32_t i, int32_t len) const
01213     {
01214         Qfloat *data;
01215         int32_t start;
01216         if((start = cache->get_data(i,&data,len)) < len)
01217         {
01218             for(int32_t j=start;j<len;j++)
01219                 data[j] = (Qfloat) kernel_function(i,j);
01220         }
01221         return data;
01222     }
01223 
01224     Qfloat *get_QD() const
01225     {
01226         return QD;
01227     }
01228 
01229     void swap_index(int32_t i, int32_t j) const
01230     {
01231         cache->swap_index(i,j);
01232         Kernel::swap_index(i,j);
01233         swap(QD[i],QD[j]);
01234     }
01235 
01236     ~ONE_CLASS_Q()
01237     {
01238         delete cache;
01239         delete[] QD;
01240     }
01241 private:
01242     Cache *cache;
01243     Qfloat *QD;
01244 };
01245 
01246 class SVR_Q: public Kernel
01247 { 
01248 public:
01249     SVR_Q(const svm_problem& prob, const svm_parameter& param)
01250     :Kernel(prob.l, prob.x, param)
01251     {
01252         l = prob.l;
01253         cache = new Cache(l,(int64_t)(param.cache_size*(1l<<20)));
01254         QD = new Qfloat[2*l];
01255         sign = new schar[2*l];
01256         index = new int32_t[2*l];
01257         for(int32_t k=0;k<l;k++)
01258         {
01259             sign[k] = 1;
01260             sign[k+l] = -1;
01261             index[k] = k;
01262             index[k+l] = k;
01263             QD[k]= (Qfloat)kernel_function(k,k);
01264             QD[k+l]=QD[k];
01265         }
01266         buffer[0] = new Qfloat[2*l];
01267         buffer[1] = new Qfloat[2*l];
01268         next_buffer = 0;
01269     }
01270 
01271     void swap_index(int32_t i, int32_t j) const
01272     {
01273         swap(sign[i],sign[j]);
01274         swap(index[i],index[j]);
01275         swap(QD[i],QD[j]);
01276     }
01277     
01278     Qfloat *get_Q(int32_t i, int32_t len) const
01279     {
01280         Qfloat *data;
01281         int32_t real_i = index[i];
01282         if(cache->get_data(real_i,&data,l) < l)
01283         {
01284             for(int32_t j=0;j<l;j++)
01285                 data[j] = (Qfloat)kernel_function(real_i,j);
01286         }
01287 
01288         // reorder and copy
01289         Qfloat *buf = buffer[next_buffer];
01290         next_buffer = 1 - next_buffer;
01291         schar si = sign[i];
01292         for(int32_t j=0;j<len;j++)
01293             buf[j] = si * sign[j] * data[index[j]];
01294         return buf;
01295     }
01296 
01297     Qfloat *get_QD() const
01298     {
01299         return QD;
01300     }
01301 
01302     ~SVR_Q()
01303     {
01304         delete cache;
01305         delete[] sign;
01306         delete[] index;
01307         delete[] buffer[0];
01308         delete[] buffer[1];
01309         delete[] QD;
01310     }
01311 
01312 private:
01313     int32_t l;
01314     Cache *cache;
01315     schar *sign;
01316     int32_t *index;
01317     mutable int32_t next_buffer;
01318     Qfloat *buffer[2];
01319     Qfloat *QD;
01320 };
01321 
01322 //
01323 // construct and solve various formulations
01324 //
01325 static void solve_c_svc(
01326     const svm_problem *prob, const svm_parameter* param,
01327     float64_t *alpha, Solver::SolutionInfo* si, float64_t Cp, float64_t Cn)
01328 {
01329     int32_t l = prob->l;
01330     float64_t *minus_ones = new float64_t[l];
01331     schar *y = new schar[l];
01332 
01333     int32_t i;
01334 
01335     for(i=0;i<l;i++)
01336     {
01337         alpha[i] = 0;
01338         minus_ones[i] = -1;
01339         if(prob->y[i] > 0) y[i] = +1; else y[i]=-1;
01340     }
01341 
01342     Solver s;
01343     s.Solve(l, SVC_Q(*prob,*param,y), minus_ones, y,
01344         alpha, Cp, Cn, param->eps, si, param->shrinking);
01345 
01346     float64_t sum_alpha=0;
01347     for(i=0;i<l;i++)
01348         sum_alpha += alpha[i];
01349 
01350     if (Cp==Cn)
01351         SG_SINFO("nu = %f\n", sum_alpha/(param->C*prob->l));
01352 
01353     for(i=0;i<l;i++)
01354         alpha[i] *= y[i];
01355 
01356     delete[] minus_ones;
01357     delete[] y;
01358 }
01359 
01360 static void solve_nu_svc(
01361     const svm_problem *prob, const svm_parameter *param,
01362     float64_t *alpha, Solver::SolutionInfo* si)
01363 {
01364     int32_t i;
01365     int32_t l = prob->l;
01366     float64_t nu = param->nu;
01367 
01368     schar *y = new schar[l];
01369 
01370     for(i=0;i<l;i++)
01371         if(prob->y[i]>0)
01372             y[i] = +1;
01373         else
01374             y[i] = -1;
01375 
01376     float64_t sum_pos = nu*l/2;
01377     float64_t sum_neg = nu*l/2;
01378 
01379     for(i=0;i<l;i++)
01380         if(y[i] == +1)
01381         {
01382             alpha[i] = min(1.0,sum_pos);
01383             sum_pos -= alpha[i];
01384         }
01385         else
01386         {
01387             alpha[i] = min(1.0,sum_neg);
01388             sum_neg -= alpha[i];
01389         }
01390 
01391     float64_t *zeros = new float64_t[l];
01392 
01393     for(i=0;i<l;i++)
01394         zeros[i] = 0;
01395 
01396     Solver_NU s;
01397     s.Solve(l, SVC_Q(*prob,*param,y), zeros, y,
01398         alpha, 1.0, 1.0, param->eps, si,  param->shrinking);
01399     float64_t r = si->r;
01400 
01401     SG_SINFO("C = %f\n",1/r);
01402 
01403     for(i=0;i<l;i++)
01404         alpha[i] *= y[i]/r;
01405 
01406     si->rho /= r;
01407     si->obj /= (r*r);
01408     si->upper_bound_p = 1/r;
01409     si->upper_bound_n = 1/r;
01410 
01411     delete[] y;
01412     delete[] zeros;
01413 }
01414 
01415 static void solve_one_class(
01416     const svm_problem *prob, const svm_parameter *param,
01417     float64_t *alpha, Solver::SolutionInfo* si)
01418 {
01419     int32_t l = prob->l;
01420     float64_t *zeros = new float64_t[l];
01421     schar *ones = new schar[l];
01422     int32_t i;
01423 
01424     int32_t n = (int32_t)(param->nu*prob->l);   // # of alpha's at upper bound
01425 
01426     for(i=0;i<n;i++)
01427         alpha[i] = 1;
01428     if(n<prob->l)
01429         alpha[n] = param->nu * prob->l - n;
01430     for(i=n+1;i<l;i++)
01431         alpha[i] = 0;
01432 
01433     for(i=0;i<l;i++)
01434     {
01435         zeros[i] = 0;
01436         ones[i] = 1;
01437     }
01438 
01439     Solver s;
01440     s.Solve(l, ONE_CLASS_Q(*prob,*param), zeros, ones,
01441         alpha, 1.0, 1.0, param->eps, si, param->shrinking);
01442 
01443     delete[] zeros;
01444     delete[] ones;
01445 }
01446 
01447 static void solve_epsilon_svr(
01448     const svm_problem *prob, const svm_parameter *param,
01449     float64_t *alpha, Solver::SolutionInfo* si)
01450 {
01451     int32_t l = prob->l;
01452     float64_t *alpha2 = new float64_t[2*l];
01453     float64_t *linear_term = new float64_t[2*l];
01454     schar *y = new schar[2*l];
01455     int32_t i;
01456 
01457     for(i=0;i<l;i++)
01458     {
01459         alpha2[i] = 0;
01460         linear_term[i] = param->p - prob->y[i];
01461         y[i] = 1;
01462 
01463         alpha2[i+l] = 0;
01464         linear_term[i+l] = param->p + prob->y[i];
01465         y[i+l] = -1;
01466     }
01467 
01468     Solver s;
01469     s.Solve(2*l, SVR_Q(*prob,*param), linear_term, y,
01470         alpha2, param->C, param->C, param->eps, si, param->shrinking);
01471 
01472     float64_t sum_alpha = 0;
01473     for(i=0;i<l;i++)
01474     {
01475         alpha[i] = alpha2[i] - alpha2[i+l];
01476         sum_alpha += fabs(alpha[i]);
01477     }
01478     SG_SINFO("nu = %f\n",sum_alpha/(param->C*l));
01479 
01480     delete[] alpha2;
01481     delete[] linear_term;
01482     delete[] y;
01483 }
01484 
01485 static void solve_nu_svr(
01486     const svm_problem *prob, const svm_parameter *param,
01487     float64_t *alpha, Solver::SolutionInfo* si)
01488 {
01489     int32_t l = prob->l;
01490     float64_t C = param->C;
01491     float64_t *alpha2 = new float64_t[2*l];
01492     float64_t *linear_term = new float64_t[2*l];
01493     schar *y = new schar[2*l];
01494     int32_t i;
01495 
01496     float64_t sum = C * param->nu * l / 2;
01497     for(i=0;i<l;i++)
01498     {
01499         alpha2[i] = alpha2[i+l] = min(sum,C);
01500         sum -= alpha2[i];
01501 
01502         linear_term[i] = - prob->y[i];
01503         y[i] = 1;
01504 
01505         linear_term[i+l] = prob->y[i];
01506         y[i+l] = -1;
01507     }
01508 
01509     Solver_NU s;
01510     s.Solve(2*l, SVR_Q(*prob,*param), linear_term, y,
01511         alpha2, C, C, param->eps, si, param->shrinking);
01512 
01513     SG_SINFO("epsilon = %f\n",-si->r);
01514 
01515     for(i=0;i<l;i++)
01516         alpha[i] = alpha2[i] - alpha2[i+l];
01517 
01518     delete[] alpha2;
01519     delete[] linear_term;
01520     delete[] y;
01521 }
01522 
01523 //
01524 // decision_function
01525 //
01526 struct decision_function
01527 {
01528     float64_t *alpha;
01529     float64_t rho;  
01530     float64_t objective;
01531 };
01532 
01533 decision_function svm_train_one(
01534     const svm_problem *prob, const svm_parameter *param,
01535     float64_t Cp, float64_t Cn)
01536 {
01537     float64_t *alpha = Malloc(float64_t, prob->l);
01538     Solver::SolutionInfo si;
01539     switch(param->svm_type)
01540     {
01541         case C_SVC:
01542             solve_c_svc(prob,param,alpha,&si,Cp,Cn);
01543             break;
01544         case NU_SVC:
01545             solve_nu_svc(prob,param,alpha,&si);
01546             break;
01547         case ONE_CLASS:
01548             solve_one_class(prob,param,alpha,&si);
01549             break;
01550         case EPSILON_SVR:
01551             solve_epsilon_svr(prob,param,alpha,&si);
01552             break;
01553         case NU_SVR:
01554             solve_nu_svr(prob,param,alpha,&si);
01555             break;
01556     }
01557 
01558     SG_SINFO("obj = %.16f, rho = %.16f\n",si.obj,si.rho);
01559 
01560     // output SVs
01561     if (param->svm_type != ONE_CLASS)
01562     {
01563         int32_t nSV = 0;
01564         int32_t nBSV = 0;
01565         for(int32_t i=0;i<prob->l;i++)
01566         {
01567             if(fabs(alpha[i]) > 0)
01568             {
01569                 ++nSV;
01570                 if(prob->y[i] > 0)
01571                 {
01572                     if(fabs(alpha[i]) >= si.upper_bound_p)
01573                         ++nBSV;
01574                 }
01575                 else
01576                 {
01577                     if(fabs(alpha[i]) >= si.upper_bound_n)
01578                         ++nBSV;
01579                 }
01580             }
01581         }
01582         SG_SINFO("nSV = %d, nBSV = %d\n",nSV,nBSV);
01583     }
01584 
01585     decision_function f;
01586     f.alpha = alpha;
01587     f.rho = si.rho;
01588     f.objective=si.obj;
01589     return f;
01590 }
01591 
01592 //
01593 // label: label name, start: begin of each class, count: #data of classes, perm: indices to the original data
01594 // perm, length l, must be allocated before calling this subroutine
01595 void svm_group_classes(
01596     const svm_problem *prob, int32_t *nr_class_ret, int32_t **label_ret,
01597     int32_t **start_ret, int32_t **count_ret, int32_t *perm)
01598 {
01599     int32_t l = prob->l;
01600     int32_t max_nr_class = 16;
01601     int32_t nr_class = 0;
01602     int32_t *label = Malloc(int32_t, max_nr_class);
01603     int32_t *count = Malloc(int32_t, max_nr_class);
01604     int32_t *data_label = Malloc(int32_t, l);
01605     int32_t i;
01606 
01607     for(i=0;i<l;i++)
01608     {
01609         int32_t this_label=(int32_t) prob->y[i];
01610         int32_t j;
01611         for(j=0;j<nr_class;j++)
01612         {
01613             if(this_label == label[j])
01614             {
01615                 ++count[j];
01616                 break;
01617             }
01618         }
01619         data_label[i] = j;
01620         if(j == nr_class)
01621         {
01622             if(nr_class == max_nr_class)
01623             {
01624                 max_nr_class *= 2;
01625                 label=(int32_t *) realloc(label,max_nr_class*sizeof(int32_t));
01626                 count=(int32_t *) realloc(count,max_nr_class*sizeof(int32_t));
01627             }
01628             label[nr_class] = this_label;
01629             count[nr_class] = 1;
01630             ++nr_class;
01631         }
01632     }
01633 
01634     int32_t *start = Malloc(int32_t, nr_class);
01635     start[0] = 0;
01636     for(i=1;i<nr_class;i++)
01637         start[i] = start[i-1]+count[i-1];
01638     for(i=0;i<l;i++)
01639     {
01640         perm[start[data_label[i]]] = i;
01641         ++start[data_label[i]];
01642     }
01643     start[0] = 0;
01644     for(i=1;i<nr_class;i++)
01645         start[i] = start[i-1]+count[i-1];
01646 
01647     *nr_class_ret = nr_class;
01648     *label_ret = label;
01649     *start_ret = start;
01650     *count_ret = count;
01651     free(data_label);
01652 }
01653 
01654 //
01655 // Interface functions
01656 //
01657 svm_model *svm_train(const svm_problem *prob, const svm_parameter *param)
01658 {
01659     svm_model *model = Malloc(svm_model,1);
01660     model->param = *param;
01661     model->free_sv = 0; // XXX
01662 
01663     if(param->svm_type == ONE_CLASS ||
01664        param->svm_type == EPSILON_SVR ||
01665        param->svm_type == NU_SVR)
01666     {
01667         SG_SINFO("training one class svm or doing epsilon sv regression\n");
01668 
01669         // regression or one-class-svm
01670         model->nr_class = 2;
01671         model->label = NULL;
01672         model->nSV = NULL;
01673         model->sv_coef = Malloc(float64_t *,1);
01674         decision_function f = svm_train_one(prob,param,0,0);
01675         model->rho = Malloc(float64_t, 1);
01676         model->rho[0] = f.rho;
01677         model->objective = f.objective;
01678 
01679         int32_t nSV = 0;
01680         int32_t i;
01681         for(i=0;i<prob->l;i++)
01682             if(fabs(f.alpha[i]) > 0) ++nSV;
01683         model->l = nSV;
01684         model->SV = Malloc(svm_node *,nSV);
01685         model->sv_coef[0] = Malloc(float64_t, nSV);
01686         int32_t j = 0;
01687         for(i=0;i<prob->l;i++)
01688             if(fabs(f.alpha[i]) > 0)
01689             {
01690                 model->SV[j] = prob->x[i];
01691                 model->sv_coef[0][j] = f.alpha[i];
01692                 ++j;
01693             }       
01694 
01695         free(f.alpha);
01696     }
01697     else
01698     {
01699         // classification
01700         int32_t l = prob->l;
01701         int32_t nr_class;
01702         int32_t *label = NULL;
01703         int32_t *start = NULL;
01704         int32_t *count = NULL;
01705         int32_t *perm = Malloc(int32_t, l);
01706 
01707         // group training data of the same class
01708         svm_group_classes(prob,&nr_class,&label,&start,&count,perm);
01709         svm_node **x = Malloc(svm_node *,l);
01710         int32_t i;
01711         for(i=0;i<l;i++)
01712             x[i] = prob->x[perm[i]];
01713 
01714         // calculate weighted C
01715 
01716         float64_t *weighted_C = Malloc(float64_t,  nr_class);
01717         for(i=0;i<nr_class;i++)
01718             weighted_C[i] = param->C;
01719         for(i=0;i<param->nr_weight;i++)
01720         {
01721             int32_t j;
01722             for(j=0;j<nr_class;j++)
01723                 if(param->weight_label[i] == label[j])
01724                     break;
01725             if(j == nr_class)
01726                 fprintf(stderr,"warning: class label %d specified in weight is not found\n", param->weight_label[i]);
01727             else
01728                 weighted_C[j] *= param->weight[i];
01729         }
01730 
01731         // train k*(k-1)/2 models
01732         
01733         bool *nonzero = Malloc(bool,l);
01734         for(i=0;i<l;i++)
01735             nonzero[i] = false;
01736         decision_function *f = Malloc(decision_function,nr_class*(nr_class-1)/2);
01737 
01738         int32_t p = 0;
01739         for(i=0;i<nr_class;i++)
01740             for(int32_t j=i+1;j<nr_class;j++)
01741             {
01742                 svm_problem sub_prob;
01743                 int32_t si = start[i], sj = start[j];
01744                 int32_t ci = count[i], cj = count[j];
01745                 sub_prob.l = ci+cj;
01746                 sub_prob.x = Malloc(svm_node *,sub_prob.l);
01747                 sub_prob.y = Malloc(float64_t, sub_prob.l+1); //dirty hack to surpress valgrind err
01748 
01749                 int32_t k;
01750                 for(k=0;k<ci;k++)
01751                 {
01752                     sub_prob.x[k] = x[si+k];
01753                     sub_prob.y[k] = +1;
01754                 }
01755                 for(k=0;k<cj;k++)
01756                 {
01757                     sub_prob.x[ci+k] = x[sj+k];
01758                     sub_prob.y[ci+k] = -1;
01759                 }
01760                 sub_prob.y[sub_prob.l]=-1; //dirty hack to surpress valgrind err
01761                 
01762                 f[p] = svm_train_one(&sub_prob,param,weighted_C[i],weighted_C[j]);
01763                 for(k=0;k<ci;k++)
01764                     if(!nonzero[si+k] && fabs(f[p].alpha[k]) > 0)
01765                         nonzero[si+k] = true;
01766                 for(k=0;k<cj;k++)
01767                     if(!nonzero[sj+k] && fabs(f[p].alpha[ci+k]) > 0)
01768                         nonzero[sj+k] = true;
01769                 free(sub_prob.x);
01770                 free(sub_prob.y);
01771                 ++p;
01772             }
01773 
01774         // build output
01775 
01776         model->objective = f[0].objective;
01777         model->nr_class = nr_class;
01778         
01779         model->label = Malloc(int32_t, nr_class);
01780         for(i=0;i<nr_class;i++)
01781             model->label[i] = label[i];
01782         
01783         model->rho = Malloc(float64_t, nr_class*(nr_class-1)/2);
01784         for(i=0;i<nr_class*(nr_class-1)/2;i++)
01785             model->rho[i] = f[i].rho;
01786 
01787         int32_t total_sv = 0;
01788         int32_t *nz_count = Malloc(int32_t, nr_class);
01789         model->nSV = Malloc(int32_t, nr_class);
01790         for(i=0;i<nr_class;i++)
01791         {
01792             int32_t nSV = 0;
01793             for(int32_t j=0;j<count[i];j++)
01794                 if(nonzero[start[i]+j])
01795                 {   
01796                     ++nSV;
01797                     ++total_sv;
01798                 }
01799             model->nSV[i] = nSV;
01800             nz_count[i] = nSV;
01801         }
01802         
01803         SG_SINFO("Total nSV = %d\n",total_sv);
01804 
01805         model->l = total_sv;
01806         model->SV = Malloc(svm_node *,total_sv);
01807         p = 0;
01808         for(i=0;i<l;i++)
01809             if(nonzero[i]) model->SV[p++] = x[i];
01810 
01811         int32_t *nz_start = Malloc(int32_t, nr_class);
01812         nz_start[0] = 0;
01813         for(i=1;i<nr_class;i++)
01814             nz_start[i] = nz_start[i-1]+nz_count[i-1];
01815 
01816         model->sv_coef = Malloc(float64_t *,nr_class-1);
01817         for(i=0;i<nr_class-1;i++)
01818             model->sv_coef[i] = Malloc(float64_t, total_sv);
01819 
01820         p = 0;
01821         for(i=0;i<nr_class;i++)
01822             for(int32_t j=i+1;j<nr_class;j++)
01823             {
01824                 // classifier (i,j): coefficients with
01825                 // i are in sv_coef[j-1][nz_start[i]...],
01826                 // j are in sv_coef[i][nz_start[j]...]
01827 
01828                 int32_t si = start[i];
01829                 int32_t sj = start[j];
01830                 int32_t ci = count[i];
01831                 int32_t cj = count[j];
01832 
01833                 int32_t q = nz_start[i];
01834                 int32_t k;
01835                 for(k=0;k<ci;k++)
01836                     if(nonzero[si+k])
01837                         model->sv_coef[j-1][q++] = f[p].alpha[k];
01838                 q = nz_start[j];
01839                 for(k=0;k<cj;k++)
01840                     if(nonzero[sj+k])
01841                         model->sv_coef[i][q++] = f[p].alpha[ci+k];
01842                 ++p;
01843             }
01844         
01845         free(label);
01846         free(count);
01847         free(perm);
01848         free(start);
01849         free(x);
01850         free(weighted_C);
01851         free(nonzero);
01852         for(i=0;i<nr_class*(nr_class-1)/2;i++)
01853             free(f[i].alpha);
01854         free(f);
01855         free(nz_count);
01856         free(nz_start);
01857     }
01858     return model;
01859 }
01860 
01861 const char *svm_type_table[] =
01862 {
01863     "c_svc","nu_svc","one_class","epsilon_svr","nu_svr",NULL
01864 };
01865 
01866 const char *kernel_type_table[]=
01867 {
01868     "linear","polynomial","rbf","sigmoid","precomputed",NULL
01869 };
01870 
01871 void svm_destroy_model(svm_model* model)
01872 {
01873     if(model->free_sv && model->l > 0)
01874         free((void *)(model->SV[0]));
01875     for(int32_t i=0;i<model->nr_class-1;i++)
01876         free(model->sv_coef[i]);
01877     free(model->SV);
01878     free(model->sv_coef);
01879     free(model->rho);
01880     free(model->label);
01881     free(model->nSV);
01882     free(model);
01883 }
01884 
01885 void svm_destroy_param(svm_parameter* param)
01886 {
01887     free(param->weight_label);
01888     free(param->weight);
01889 }
01890 
01891 const char *svm_check_parameter(
01892     const svm_problem *prob, const svm_parameter *param)
01893 {
01894     // svm_type
01895 
01896     int32_t svm_type = param->svm_type;
01897     if(svm_type != C_SVC &&
01898        svm_type != NU_SVC &&
01899        svm_type != ONE_CLASS &&
01900        svm_type != EPSILON_SVR &&
01901        svm_type != NU_SVR)
01902         return "unknown svm type";
01903 
01904     // kernel_type, degree
01905 
01906     int32_t kernel_type = param->kernel_type;
01907     if(kernel_type != LINEAR &&
01908        kernel_type != POLY &&
01909        kernel_type != RBF &&
01910        kernel_type != SIGMOID &&
01911        kernel_type != PRECOMPUTED)
01912         return "unknown kernel type";
01913 
01914     if(param->degree < 0)
01915         return "degree of polynomial kernel < 0";
01916 
01917     // cache_size,eps,C,nu,p,shrinking
01918 
01919     if(param->cache_size <= 0)
01920         return "cache_size <= 0";
01921 
01922     if(param->eps <= 0)
01923         return "eps <= 0";
01924 
01925     if(svm_type == C_SVC ||
01926        svm_type == EPSILON_SVR ||
01927        svm_type == NU_SVR)
01928         if(param->C <= 0)
01929             return "C <= 0";
01930 
01931     if(svm_type == NU_SVC ||
01932        svm_type == ONE_CLASS ||
01933        svm_type == NU_SVR)
01934         if(param->nu <= 0 || param->nu > 1)
01935             return "nu <= 0 or nu > 1";
01936 
01937     if(svm_type == EPSILON_SVR)
01938         if(param->p < 0)
01939             return "p < 0";
01940 
01941     if(param->shrinking != 0 &&
01942        param->shrinking != 1)
01943         return "shrinking != 0 and shrinking != 1";
01944 
01945 
01946     // check whether nu-svc is feasible
01947     
01948     if(svm_type == NU_SVC)
01949     {
01950         int32_t l = prob->l;
01951         int32_t max_nr_class = 16;
01952         int32_t nr_class = 0;
01953         int32_t *label = Malloc(int32_t, max_nr_class);
01954         int32_t *count = Malloc(int32_t, max_nr_class);
01955 
01956         int32_t i;
01957         for(i=0;i<l;i++)
01958         {
01959             int32_t this_label = (int32_t) prob->y[i];
01960             int32_t j;
01961             for(j=0;j<nr_class;j++)
01962                 if(this_label == label[j])
01963                 {
01964                     ++count[j];
01965                     break;
01966                 }
01967             if(j == nr_class)
01968             {
01969                 if(nr_class == max_nr_class)
01970                 {
01971                     max_nr_class *= 2;
01972                     label=(int32_t *) realloc(label,
01973                         max_nr_class*sizeof(int32_t));
01974                     count=(int32_t *) realloc(count,
01975                         max_nr_class*sizeof(int32_t));
01976                 }
01977                 label[nr_class] = this_label;
01978                 count[nr_class] = 1;
01979                 ++nr_class;
01980             }
01981         }
01982     
01983         for(i=0;i<nr_class;i++)
01984         {
01985             int32_t n1 = count[i];
01986             for(int32_t j=i+1;j<nr_class;j++)
01987             {
01988                 int32_t n2 = count[j];
01989                 if(param->nu*(n1+n2)/2 > min(n1,n2))
01990                 {
01991                     free(label);
01992                     free(count);
01993                     return "specified nu is infeasible";
01994                 }
01995             }
01996         }
01997         free(label);
01998         free(count);
01999     }
02000 
02001     return NULL;
02002 }

SHOGUN Machine Learning Toolbox - Documentation