SVM_libsvm.cpp

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

SHOGUN Machine Learning Toolbox - Documentation