开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > Windows编程 > .net编程 > 查看源码
svm.cpp
资源名称：svm.rar　[点击查看]
上传用户：xgw_05
上传日期：2014-12-08
资源大小：2726k
文件大小：42k
源码类别：

.net编程

开发平台：

Java

svm.cpp：源码内容
  1. #include <math.h>
  2. #include <stdio.h>
  3. #include <stdlib.h>
  4. #include <ctype.h>
  5. #include <float.h>
  6. #include <string.h>
  7. #include <stdarg.h>
  8. #include "svm.h"
  9. typedef float Qfloat;
  10. typedef signed char schar;
  11. #ifndef min
  12. template <class T> inline T min(T x,T y) { return (x<y)?x:y; }
  13. #endif
  14. #ifndef max
  15. template <class T> inline T max(T x,T y) { return (x>y)?x:y; }
  16. #endif
  17. template <class T> inline void swap(T& x, T& y) { T t=x; x=y; y=t; }
  18. template <class S, class T> inline void clone(T*& dst, S* src, int n)
  19. {
  20. dst = new T[n];
  21. memcpy((void *)dst,(void *)src,sizeof(T)*n);
  22. }
  23. #define INF HUGE_VAL
  24. #define Malloc(type,n) (type *)malloc((n)*sizeof(type))
  25. #if 1
  26. void info(char *fmt,...)
  27. {
  28. va_list ap;
  29. va_start(ap,fmt);
  30. vprintf(fmt,ap);
  31. va_end(ap);
  32. }
  33. void info_flush()
  34. {
  35. fflush(stdout);
  36. }
  37. #else
  38. void info(char *fmt,...) {}
  39. void info_flush() {}
  40. #endif
  42. //
  43. // Kernel Cache
  44. //
  45. // l is the number of total data items
  46. // size is the cache size limit in bytes
  47. //
  48. class Cache
  49. {
  50. public:
  51. Cache(int l,int size);
  52. ~Cache();
  54. // request data [0,len)
  55. // return some position p where [p,len) need to be filled
  56. // (p >= len if nothing needs to be filled)
  57. int get_data(const int index, Qfloat **data, int len);
  58. void swap_index(int i, int j); // future_option
  59. private:
  60. int l;
  61. int size;
  62. struct head_t
  63. {
  64. head_t *prev, *next; // a cicular list
  65. Qfloat *data;
  66. int len; // data[0,len) is cached in this entry
  67. };
  69. head_t* head;
  70. head_t lru_head;
  71. void lru_delete(head_t *h);
  72. void lru_insert(head_t *h);
  73. };
  75. Cache::Cache(int l_,int size_):l(l_),size(size_)
  76. {
  77. head = (head_t *)calloc(l,sizeof(head_t)); // initialized to 0
  78. size /= sizeof(Qfloat);
  79. size -= l * sizeof(head_t) / sizeof(Qfloat);
  80. lru_head.next = lru_head.prev = &lru_head;
  81. }
  83. Cache::~Cache()
  84. {
  85. for(head_t *h = lru_head.next; h != &lru_head; h=h->next)
  86. free(h->data);
  87. free(head);
  88. }
  90. void Cache::lru_delete(head_t *h)
  91. {
  92. // delete from current location
  93. h->prev->next = h->next;
  94. h->next->prev = h->prev;
  95. }
  97. void Cache::lru_insert(head_t *h)
  98. {
  99. // insert to last position
  100. h->next = &lru_head;
  101. h->prev = lru_head.prev;
  102. h->prev->next = h;
  103. h->next->prev = h;
  104. }
  106. int Cache::get_data(const int index, Qfloat **data, int len)
  107. {
  108. head_t *h = &head[index];
  109. if(h->len) lru_delete(h);
  110. int more = len - h->len;
  112. if(more > 0)
  113. {
  114. // free old space
  115. while(size < more)
  116. {
  117. head_t *old = lru_head.next;
  118. lru_delete(old);
  119. free(old->data);
  120. size += old->len;
  121. old->data = 0;
  122. old->len = 0;
  123. }
  125. // allocate new space
  126. h->data = (Qfloat *)realloc(h->data,sizeof(Qfloat)*len);
  127. size -= more;
  128. swap(h->len,len);
  129. }
  131. lru_insert(h);
  132. *data = h->data;
  133. return len;
  134. }
  136. void Cache::swap_index(int i, int j)
  137. {
  138. if(i==j) return;
  140. if(head[i].len) lru_delete(&head[i]);
  141. if(head[j].len) lru_delete(&head[j]);
  142. swap(head[i].data,head[j].data);
  143. swap(head[i].len,head[j].len);
  144. if(head[i].len) lru_insert(&head[i]);
  145. if(head[j].len) lru_insert(&head[j]);
  147. if(i>j) swap(i,j);
  148. for(head_t *h = lru_head.next; h!=&lru_head; h=h->next)
  149. {
  150. if(h->len > i)
  151. {
  152. if(h->len > j)
  153. swap(h->data[i],h->data[j]);
  154. else
  155. {
  156. // give up
  157. lru_delete(h);
  158. free(h->data);
  159. size += h->len;
  160. h->data = 0;
  161. h->len = 0;
  162. }
  163. }
  164. }
  165. }
  167. //
  168. // Kernel evaluation
  169. //
  170. // the static method k_function is for doing single kernel evaluation
  171. // the constructor of Kernel prepares to calculate the l*l kernel matrix
  172. // the member function get_Q is for getting one column from the Q Matrix
  173. //
  174. class Kernel {
  175. public:
  176. Kernel(int l, svm_node * const * x, const svm_parameter& param);
  177. virtual ~Kernel();
  179. static double k_function(const svm_node *x, const svm_node *y,
  180.  const svm_parameter& param);
  181. virtual Qfloat *get_Q(int column, int len) const = 0;
  182. virtual void swap_index(int i, int j) const // no so const...
  183. {
  184. swap(x[i],x[j]);
  185. if(x_square) swap(x_square[i],x_square[j]);
  186. }
  187. protected:
  189. double (Kernel::*kernel_function)(int i, int j) const;
  191. private:
  192. const svm_node **x;
  193. double *x_square;
  195. // svm_parameter
  196. const int kernel_type;
  197. const double degree;
  198. const double gamma;
  199. const double coef0;
  201. static double dot(const svm_node *px, const svm_node *py);
  202. double kernel_linear(int i, int j) const
  203. {
  204. return dot(x[i],x[j]);
  205. }
  206. double kernel_poly(int i, int j) const
  207. {
  208. return pow(gamma*dot(x[i],x[j])+coef0,degree);
  209. }
  210. double kernel_rbf(int i, int j) const
  211. {
  212. return exp(-gamma*(x_square[i]+x_square[j]-2*dot(x[i],x[j])));
  213. }
  214. double kernel_sigmoid(int i, int j) const
  215. {
  216. return tanh(gamma*dot(x[i],x[j])+coef0);
  217. }
  218. };
  220. Kernel::Kernel(int l, svm_node * const * x_, const svm_parameter& param)
  221. :kernel_type(param.kernel_type), degree(param.degree),
  222.  gamma(param.gamma), coef0(param.coef0)
  223. {
  224. switch(kernel_type)
  225. {
  226. case LINEAR:
  227. kernel_function = &Kernel::kernel_linear;
  228. break;
  229. case POLY:
  230. kernel_function = &Kernel::kernel_poly;
  231. break;
  232. case RBF:
  233. kernel_function = &Kernel::kernel_rbf;
  234. break;
  235. case SIGMOID:
  236. kernel_function = &Kernel::kernel_sigmoid;
  237. break;
  238. default:
  239. fprintf(stderr,"unknown kernel function.n");
  240. exit(1);
  241. }
  243. clone(x,x_,l);
  245. if(kernel_type == RBF)
  246. {
  247. x_square = new double[l];
  248. for(int i=0;i<l;i++)
  249. x_square[i] = dot(x[i],x[i]);
  250. }
  251. else
  252. x_square = 0;
  253. }
  255. Kernel::~Kernel()
  256. {
  257. delete[] x;
  258. delete[] x_square;
  259. }
  261. double Kernel::dot(const svm_node *px, const svm_node *py)
  262. {
  263. double sum = 0;
  264. while(px->index != -1 && py->index != -1)
  265. {
  266. if(px->index == py->index)
  267. {
  268. sum += px->value * py->value;
  269. ++px;
  270. ++py;
  271. }
  272. else
  273. {
  274. if(px->index > py->index)
  275. ++py;
  276. else
  277. ++px;
  278. }
  279. }
  280. return sum;
  281. }
  283. double Kernel::k_function(const svm_node *x, const svm_node *y,
  284.   const svm_parameter& param)
  285. {
  286. switch(param.kernel_type)
  287. {
  288. case LINEAR:
  289. return dot(x,y);
  290. case POLY:
  291. return pow(param.gamma*dot(x,y)+param.coef0,param.degree);
  292. case RBF:
  293. {
  294. double sum = 0;
  295. while(x->index != -1 && y->index !=-1)
  296. {
  297. if(x->index == y->index)
  298. {
  299. double d = x->value - y->value;
  300. sum += d*d;
  301. ++x;
  302. ++y;
  303. }
  304. else
  305. {
  306. if(x->index > y->index)
  307. {
  308. sum += y->value * y->value;
  309. ++y;
  310. }
  311. else
  312. {
  313. sum += x->value * x->value;
  314. ++x;
  315. }
  316. }
  317. }
  319. while(x->index != -1)
  320. {
  321. sum += x->value * x->value;
  322. ++x;
  323. }
  325. while(y->index != -1)
  326. {
  327. sum += y->value * y->value;
  328. ++y;
  329. }
  331. return exp(-param.gamma*sum);
  332. }
  333. case SIGMOID:
  334. return tanh(param.gamma*dot(x,y)+param.coef0);
  335. default:
  336. break;
  337. }
  338. fprintf(stderr,"unknown kernel function.n");
  339. exit(1);
  340. }
  342. // Generalized SMO+SVMlight algorithm
  343. // Solves:
  344. //
  345. // min 0.5(alpha^T Q alpha) + b^T alpha
  346. //
  347. // y^T alpha = delta
  348. // y_i = +1 or -1
  349. // 0 <= alpha_i <= Cp for y_i = 1
  350. // 0 <= alpha_i <= Cn for y_i = -1
  351. //
  352. // Given:
  353. //
  354. // Q, b, y, Cp, Cn, and an initial feasible point alpha
  355. // l is the size of vectors and matrices
  356. // eps is the stopping criterion
  357. //
  358. // solution will be put in alpha, objective value will be put in obj
  359. //
  360. class Solver {
  361. public:
  362. Solver() {};
  363. virtual ~Solver() {};
  365. struct SolutionInfo {
  366. double obj;
  367. double rho;
  368. double upper_bound_p;
  369. double upper_bound_n;
  370. double r; // for Solver_NU
  371. };
  373. void Solve(int l, const Kernel& Q, const double *b_, const schar *y_,
  374.    double *alpha_, double Cp, double Cn, double eps,
  375.    SolutionInfo* si, int shrinking);
  376. protected:
  377. int active_size;
  378. schar *y;
  379. double *G; // gradient of objective function
  380. enum { LOWER_BOUND, UPPER_BOUND, FREE };
  381. char *alpha_status; // LOWER_BOUND, UPPER_BOUND, FREE
  382. double *alpha;
  383. const Kernel *Q;
  384. double eps;
  385. double Cp,Cn;
  386. double *b;
  387. int *active_set;
  388. double *G_bar; // gradient, if we treat free variables as 0
  389. int l;
  390. bool unshrinked; // XXX
  392. double get_C(int i)
  393. {
  394. return (y[i] > 0)? Cp : Cn;
  395. }
  396. void update_alpha_status(int i)
  397. {
  398. if(alpha[i] >= get_C(i))
  399. alpha_status[i] = UPPER_BOUND;
  400. else if(alpha[i] <= 0)
  401. alpha_status[i] = LOWER_BOUND;
  402. else alpha_status[i] = FREE;
  403. }
  404. bool is_upper_bound(int i) { return alpha_status[i] == UPPER_BOUND; }
  405. bool is_lower_bound(int i) { return alpha_status[i] == LOWER_BOUND; }
  406. bool is_free(int i) { return alpha_status[i] == FREE; }
  407. void swap_index(int i, int j);
  408. void reconstruct_gradient();
  409. virtual int select_working_set(int &i, int &j);
  410. virtual double calculate_rho();
  411. virtual void do_shrinking();
  412. };
  414. void Solver::swap_index(int i, int j)
  415. {
  416. Q->swap_index(i,j);
  417. swap(y[i],y[j]);
  418. swap(G[i],G[j]);
  419. swap(alpha_status[i],alpha_status[j]);
  420. swap(alpha[i],alpha[j]);
  421. swap(b[i],b[j]);
  422. swap(active_set[i],active_set[j]);
  423. swap(G_bar[i],G_bar[j]);
  424. }
  426. void Solver::reconstruct_gradient()
  427. {
  428. // reconstruct inactive elements of G from G_bar and free variables
  430. if(active_size == l) return;
  432. int i;
  433. for(i=active_size;i<l;i++)
  434. G[i] = G_bar[i] + b[i];
  436. for(i=0;i<active_size;i++)
  437. if(is_free(i))
  438. {
  439. const Qfloat *Q_i = Q->get_Q(i,l);
  440. double alpha_i = alpha[i];
  441. for(int j=active_size;j<l;j++)
  442. G[j] += alpha_i * Q_i[j];
  443. }
  444. }
  446. void Solver::Solve(int l, const Kernel& Q, const double *b_, const schar *y_,
  447.    double *alpha_, double Cp, double Cn, double eps,
  448.    SolutionInfo* si, int shrinking)
  449. {
  450. this->l = l;
  451. this->Q = &Q;
  452. clone(b, b_,l);
  453. clone(y, y_,l);
  454. clone(alpha,alpha_,l);
  455. this->Cp = Cp;
  456. this->Cn = Cn;
  457. this->eps = eps;
  458. unshrinked = false;
  460. // initialize alpha_status
  461. {
  462. alpha_status = new char[l];
  463. for(int i=0;i<l;i++)
  464. update_alpha_status(i);
  465. }
  467. // initialize active set (for shrinking)
  468. {
  469. active_set = new int[l];
  470. for(int i=0;i<l;i++)
  471. active_set[i] = i;
  472. active_size = l;
  473. }
  475. // initialize gradient
  476. {
  477. G = new double[l];
  478. G_bar = new double[l];
  479. int i;
  480. for(i=0;i<l;i++)
  481. {
  482. G[i] = b[i];
  483. G_bar[i] = 0;
  484. }
  485. for(i=0;i<l;i++)
  486. if(!is_lower_bound(i))
  487. {
  488. Qfloat *Q_i = Q.get_Q(i,l);
  489. double alpha_i = alpha[i];
  490. int j;
  491. for(j=0;j<l;j++)
  492. G[j] += alpha_i*Q_i[j];
  493. if(is_upper_bound(i))
  494. for(j=0;j<l;j++)
  495. G_bar[j] += get_C(i) * Q_i[j];
  496. }
  497. }
  499. // optimization step
  501. int iter = 0;
  502. int counter = min(l,1000)+1;
  504. while(1)
  505. {
  506. // show progress and do shrinking
  508. if(--counter == 0)
  509. {
  510. counter = min(l,1000);
  511. if(shrinking) do_shrinking();
  512. info("."); info_flush();
  513. }
  515. int i,j;
  516. if(select_working_set(i,j)!=0)
  517. {
  518. // reconstruct the whole gradient
  519. reconstruct_gradient();
  520. // reset active set size and check
  521. active_size = l;
  522. info("*"); info_flush();
  523. if(select_working_set(i,j)!=0)
  524. break;
  525. else
  526. counter = 1; // do shrinking next iteration
  527. }
  529. ++iter;
  531. // update alpha[i] and alpha[j], handle bounds carefully
  533. const Qfloat *Q_i = Q.get_Q(i,active_size);
  534. const Qfloat *Q_j = Q.get_Q(j,active_size);
  536. double C_i = get_C(i);
  537. double C_j = get_C(j);
  539. double old_alpha_i = alpha[i];
  540. double old_alpha_j = alpha[j];
  542. if(y[i]!=y[j])
  543. {
  544. double delta = (-G[i]-G[j])/(Q_i[i]+Q_j[j]+2*Q_i[j]);
  545. double diff = alpha[i] - alpha[j];
  546. alpha[i] += delta;
  547. alpha[j] += delta;
  549. if(diff > 0)
  550. {
  551. if(alpha[j] < 0)
  552. {
  553. alpha[j] = 0;
  554. alpha[i] = diff;
  555. }
  556. }
  557. else
  558. {
  559. if(alpha[i] < 0)
  560. {
  561. alpha[i] = 0;
  562. alpha[j] = -diff;
  563. }
  564. }
  565. if(diff > C_i - C_j)
  566. {
  567. if(alpha[i] > C_i)
  568. {
  569. alpha[i] = C_i;
  570. alpha[j] = C_i - diff;
  571. }
  572. }
  573. else
  574. {
  575. if(alpha[j] > C_j)
  576. {
  577. alpha[j] = C_j;
  578. alpha[i] = C_j + diff;
  579. }
  580. }
  581. }
  582. else
  583. {
  584. double delta = (G[i]-G[j])/(Q_i[i]+Q_j[j]-2*Q_i[j]);
  585. double sum = alpha[i] + alpha[j];
  586. alpha[i] -= delta;
  587. alpha[j] += delta;
  588. if(sum > C_i)
  589. {
  590. if(alpha[i] > C_i)
  591. {
  592. alpha[i] = C_i;
  593. alpha[j] = sum - C_i;
  594. }
  595. }
  596. else
  597. {
  598. if(alpha[j] < 0)
  599. {
  600. alpha[j] = 0;
  601. alpha[i] = sum;
  602. }
  603. }
  604. if(sum > C_j)
  605. {
  606. if(alpha[j] > C_j)
  607. {
  608. alpha[j] = C_j;
  609. alpha[i] = sum - C_j;
  610. }
  611. }
  612. else
  613. {
  614. if(alpha[i] < 0)
  615. {
  616. alpha[i] = 0;
  617. alpha[j] = sum;
  618. }
  619. }
  620. }
  622. // update G
  624. double delta_alpha_i = alpha[i] - old_alpha_i;
  625. double delta_alpha_j = alpha[j] - old_alpha_j;
  627. for(int k=0;k<active_size;k++)
  628. {
  629. G[k] += Q_i[k]*delta_alpha_i + Q_j[k]*delta_alpha_j;
  630. }
  632. // update alpha_status and G_bar
  634. {
  635. bool ui = is_upper_bound(i);
  636. bool uj = is_upper_bound(j);
  637. update_alpha_status(i);
  638. update_alpha_status(j);
  639. int k;
  640. if(ui != is_upper_bound(i))
  641. {
  642. Q_i = Q.get_Q(i,l);
  643. if(ui)
  644. for(k=0;k<l;k++)
  645. G_bar[k] -= C_i * Q_i[k];
  646. else
  647. for(k=0;k<l;k++)
  648. G_bar[k] += C_i * Q_i[k];
  649. }
  651. if(uj != is_upper_bound(j))
  652. {
  653. Q_j = Q.get_Q(j,l);
  654. if(uj)
  655. for(k=0;k<l;k++)
  656. G_bar[k] -= C_j * Q_j[k];
  657. else
  658. for(k=0;k<l;k++)
  659. G_bar[k] += C_j * Q_j[k];
  660. }
  661. }
  662. }
  664. // calculate rho
  666. si->rho = calculate_rho();
  668. // calculate objective value
  669. {
  670. double v = 0;
  671. int i;
  672. for(i=0;i<l;i++)
  673. v += alpha[i] * (G[i] + b[i]);
  675. si->obj = v/2;
  676. }
  678. // put back the solution
  679. {
  680. for(int i=0;i<l;i++)
  681. alpha_[active_set[i]] = alpha[i];
  682. }
  684. // juggle everything back
  685. /*{
  686. for(int i=0;i<l;i++)
  687. while(active_set[i] != i)
  688. swap_index(i,active_set[i]);
  689. // or Q.swap_index(i,active_set[i]);
  690. }*/
  692. si->upper_bound_p = Cp;
  693. si->upper_bound_n = Cn;
  695. info("noptimization finished, #iter = %dn",iter);
  697. delete[] b;
  698. delete[] y;
  699. delete[] alpha;
  700. delete[] alpha_status;
  701. delete[] active_set;
  702. delete[] G;
  703. delete[] G_bar;
  704. }
  706. // return 1 if already optimal, return 0 otherwise
  707. int Solver::select_working_set(int &out_i, int &out_j)
  708. {
  709. // return i,j which maximize -grad(f)^T d , under constraint
  710. // if alpha_i == C, d != +1
  711. // if alpha_i == 0, d != -1
  713. double Gmax1 = -INF; // max { -grad(f)_i * d | y_i*d = +1 }
  714. int Gmax1_idx = -1;
  716. double Gmax2 = -INF; // max { -grad(f)_i * d | y_i*d = -1 }
  717. int Gmax2_idx = -1;
  719. for(int i=0;i<active_size;i++)
  720. {
  721. if(y[i]==+1) // y = +1
  722. {
  723. if(!is_upper_bound(i)) // d = +1
  724. {
  725. if(-G[i] > Gmax1)
  726. {
  727. Gmax1 = -G[i];
  728. Gmax1_idx = i;
  729. }
  730. }
  731. if(!is_lower_bound(i)) // d = -1
  732. {
  733. if(G[i] > Gmax2)
  734. {
  735. Gmax2 = G[i];
  736. Gmax2_idx = i;
  737. }
  738. }
  739. }
  740. else // y = -1
  741. {
  742. if(!is_upper_bound(i)) // d = +1
  743. {
  744. if(-G[i] > Gmax2)
  745. {
  746. Gmax2 = -G[i];
  747. Gmax2_idx = i;
  748. }
  749. }
  750. if(!is_lower_bound(i)) // d = -1
  751. {
  752. if(G[i] > Gmax1)
  753. {
  754. Gmax1 = G[i];
  755. Gmax1_idx = i;
  756. }
  757. }
  758. }
  759. }
  761. if(Gmax1+Gmax2 < eps)
  762.   return 1;
  764. out_i = Gmax1_idx;
  765. out_j = Gmax2_idx;
  766. return 0;
  767. }
  769. void Solver::do_shrinking()
  770. {
  771. int i,j,k;
  772. if(select_working_set(i,j)!=0) return;
  773. double Gm1 = -y[j]*G[j];
  774. double Gm2 = y[i]*G[i];
  776. // shrink
  778. for(k=0;k<active_size;k++)
  779. {
  780. if(is_lower_bound(k))
  781. {
  782. if(y[k]==+1)
  783. {
  784. if(-G[k] >= Gm1) continue;
  785. }
  786. else if(-G[k] >= Gm2) continue;
  787. }
  788. else if(is_upper_bound(k))
  789. {
  790. if(y[k]==+1)
  791. {
  792. if(G[k] >= Gm2) continue;
  793. }
  794. else if(G[k] >= Gm1) continue;
  795. }
  796. else continue;
  798. --active_size;
  799. swap_index(k,active_size);
  800. --k; // look at the newcomer
  801. }
  803. // unshrink, check all variables again before final iterations
  805. if(unshrinked || -(Gm1 + Gm2) > eps*10) return;
  807. unshrinked = true;
  808. reconstruct_gradient();
  810. for(k=l-1;k>=active_size;k--)
  811. {
  812. if(is_lower_bound(k))
  813. {
  814. if(y[k]==+1)
  815. {
  816. if(-G[k] < Gm1) continue;
  817. }
  818. else if(-G[k] < Gm2) continue;
  819. }
  820. else if(is_upper_bound(k))
  821. {
  822. if(y[k]==+1)
  823. {
  824. if(G[k] < Gm2) continue;
  825. }
  826. else if(G[k] < Gm1) continue;
  827. }
  828. else continue;
  830. swap_index(k,active_size);
  831. active_size++;
  832. ++k; // look at the newcomer
  833. }
  834. }
  836. double Solver::calculate_rho()
  837. {
  838. double r;
  839. int nr_free = 0;
  840. double ub = INF, lb = -INF, sum_free = 0;
  841. for(int i=0;i<active_size;i++)
  842. {
  843. double yG = y[i]*G[i];
  845. if(is_lower_bound(i))
  846. {
  847. if(y[i] > 0)
  848. ub = min(ub,yG);
  849. else
  850. lb = max(lb,yG);
  851. }
  852. else if(is_upper_bound(i))
  853. {
  854. if(y[i] < 0)
  855. ub = min(ub,yG);
  856. else
  857. lb = max(lb,yG);
  858. }
  859. else
  860. {
  861. ++nr_free;
  862. sum_free += yG;
  863. }
  864. }
  866. if(nr_free>0)
  867. r = sum_free/nr_free;
  868. else
  869. r = (ub+lb)/2;
  871. return r;
  872. }
  874. //
  875. // Solver for nu-svm classification and regression
  876. //
  877. // additional constraint: e^T alpha = constant
  878. //
  879. class Solver_NU : public Solver
  880. {
  881. public:
  882. Solver_NU() {}
  883. void Solve(int l, const Kernel& Q, const double *b, const schar *y,
  884.    double *alpha, double Cp, double Cn, double eps,
  885.    SolutionInfo* si, int shrinking)
  886. {
  887. this->si = si;
  888. Solver::Solve(l,Q,b,y,alpha,Cp,Cn,eps,si,shrinking);
  889. }
  890. private:
  891. SolutionInfo *si;
  892. int select_working_set(int &i, int &j);
  893. double calculate_rho();
  894. void do_shrinking();
  895. };
  897. int Solver_NU::select_working_set(int &out_i, int &out_j)
  898. {
  899. // return i,j which maximize -grad(f)^T d , under constraint
  900. // if alpha_i == C, d != +1
  901. // if alpha_i == 0, d != -1
  903. double Gmax1 = -INF; // max { -grad(f)_i * d | y_i = +1, d = +1 }
  904. int Gmax1_idx = -1;
  906. double Gmax2 = -INF; // max { -grad(f)_i * d | y_i = +1, d = -1 }
  907. int Gmax2_idx = -1;
  909. double Gmax3 = -INF; // max { -grad(f)_i * d | y_i = -1, d = +1 }
  910. int Gmax3_idx = -1;
  912. double Gmax4 = -INF; // max { -grad(f)_i * d | y_i = -1, d = -1 }
  913. int Gmax4_idx = -1;
  915. for(int i=0;i<active_size;i++)
  916. {
  917. if(y[i]==+1) // y == +1
  918. {
  919. if(!is_upper_bound(i)) // d = +1
  920. {
  921. if(-G[i] > Gmax1)
  922. {
  923. Gmax1 = -G[i];
  924. Gmax1_idx = i;
  925. }
  926. }
  927. if(!is_lower_bound(i)) // d = -1
  928. {
  929. if(G[i] > Gmax2)
  930. {
  931. Gmax2 = G[i];
  932. Gmax2_idx = i;
  933. }
  934. }
  935. }
  936. else // y == -1
  937. {
  938. if(!is_upper_bound(i)) // d = +1
  939. {
  940. if(-G[i] > Gmax3)
  941. {
  942. Gmax3 = -G[i];
  943. Gmax3_idx = i;
  944. }
  945. }
  946. if(!is_lower_bound(i)) // d = -1
  947. {
  948. if(G[i] > Gmax4)
  949. {
  950. Gmax4 = G[i];
  951. Gmax4_idx = i;
  952. }
  953. }
  954. }
  955. }
  957. if(max(Gmax1+Gmax2,Gmax3+Gmax4) < eps)
  958.   return 1;
  960. if(Gmax1+Gmax2 > Gmax3+Gmax4)
  961. {
  962. out_i = Gmax1_idx;
  963. out_j = Gmax2_idx;
  964. }
  965. else
  966. {
  967. out_i = Gmax3_idx;
  968. out_j = Gmax4_idx;
  969. }
  970. return 0;
  971. }
  973. void Solver_NU::do_shrinking()
  974. {
  975. double Gmax1 = -INF; // max { -grad(f)_i * d | y_i = +1, d = +1 }
  976. double Gmax2 = -INF; // max { -grad(f)_i * d | y_i = +1, d = -1 }
  977. double Gmax3 = -INF; // max { -grad(f)_i * d | y_i = -1, d = +1 }
  978. double Gmax4 = -INF; // max { -grad(f)_i * d | y_i = -1, d = -1 }
  980. int k;
  981. for(k=0;k<active_size;k++)
  982. {
  983. if(!is_upper_bound(k))
  984. {
  985. if(y[k]==+1)
  986. {
  987. if(-G[k] > Gmax1) Gmax1 = -G[k];
  988. }
  989. else if(-G[k] > Gmax3) Gmax3 = -G[k];
  990. }
  991. if(!is_lower_bound(k))
  992. {
  993. if(y[k]==+1)
  994. {
  995. if(G[k] > Gmax2) Gmax2 = G[k];
  996. }
  997. else if(G[k] > Gmax4) Gmax4 = G[k];
  998. }
  999. }
  1001. double Gm1 = -Gmax2;
  1002. double Gm2 = -Gmax1;
  1003. double Gm3 = -Gmax4;
  1004. double Gm4 = -Gmax3;
  1006. for(k=0;k<active_size;k++)
  1007. {
  1008. if(is_lower_bound(k))
  1009. {
  1010. if(y[k]==+1)
  1011. {
  1012. if(-G[k] >= Gm1) continue;
  1013. }
  1014. else if(-G[k] >= Gm3) continue;
  1015. }
  1016. else if(is_upper_bound(k))
  1017. {
  1018. if(y[k]==+1)
  1019. {
  1020. if(G[k] >= Gm2) continue;
  1021. }
  1022. else if(G[k] >= Gm4) continue;
  1023. }
  1024. else continue;
  1026. --active_size;
  1027. swap_index(k,active_size);
  1028. --k; // look at the newcomer
  1029. }
  1031. // unshrink, check all variables again before final iterations
  1033. if(unshrinked || max(-(Gm1+Gm2),-(Gm3+Gm4)) > eps*10) return;
  1035. unshrinked = true;
  1036. reconstruct_gradient();
  1038. for(k=l-1;k>=active_size;k--)
  1039. {
  1040. if(is_lower_bound(k))
  1041. {
  1042. if(y[k]==+1)
  1043. {
  1044. if(-G[k] < Gm1) continue;
  1045. }
  1046. else if(-G[k] < Gm3) continue;
  1047. }
  1048. else if(is_upper_bound(k))
  1049. {
  1050. if(y[k]==+1)
  1051. {
  1052. if(G[k] < Gm2) continue;
  1053. }
  1054. else if(G[k] < Gm4) continue;
  1055. }
  1056. else continue;
  1058. swap_index(k,active_size);
  1059. active_size++;
  1060. ++k; // look at the newcomer
  1061. }
  1062. }
  1064. double Solver_NU::calculate_rho()
  1065. {
  1066. int nr_free1 = 0,nr_free2 = 0;
  1067. double ub1 = INF, ub2 = INF;
  1068. double lb1 = -INF, lb2 = -INF;
  1069. double sum_free1 = 0, sum_free2 = 0;
  1071. for(int i=0;i<active_size;i++)
  1072. {
  1073. if(y[i]==+1)
  1074. {
  1075. if(is_lower_bound(i))
  1076. ub1 = min(ub1,G[i]);
  1077. else if(is_upper_bound(i))
  1078. lb1 = max(lb1,G[i]);
  1079. else
  1080. {
  1081. ++nr_free1;
  1082. sum_free1 += G[i];
  1083. }
  1084. }
  1085. else
  1086. {
  1087. if(is_lower_bound(i))
  1088. ub2 = min(ub2,G[i]);
  1089. else if(is_upper_bound(i))
  1090. lb2 = max(lb2,G[i]);
  1091. else
  1092. {
  1093. ++nr_free2;
  1094. sum_free2 += G[i];
  1095. }
  1096. }
  1097. }
  1099. double r1,r2;
  1100. if(nr_free1 > 0)
  1101. r1 = sum_free1/nr_free1;
  1102. else
  1103. r1 = (ub1+lb1)/2;
  1105. if(nr_free2 > 0)
  1106. r2 = sum_free2/nr_free2;
  1107. else
  1108. r2 = (ub2+lb2)/2;
  1110. si->r = (r1+r2)/2;
  1111. return (r1-r2)/2;
  1112. }
  1114. //
  1115. // Q matrices for various formulations
  1116. //
  1117. class SVC_Q: public Kernel
  1118. { 
  1119. public:
  1120. SVC_Q(const svm_problem& prob, const svm_parameter& param, const schar *y_)
  1121. :Kernel(prob.l, prob.x, param)
  1122. {
  1123. clone(y,y_,prob.l);
  1124. cache = new Cache(prob.l,(int)(param.cache_size*(1<<20)));
  1125. }
  1127. Qfloat *get_Q(int i, int len) const
  1128. {
  1129. Qfloat *data;
  1130. int start;
  1131. if((start = cache->get_data(i,&data,len)) < len)
  1132. {
  1133. for(int j=start;j<len;j++)
  1134. data[j] = (Qfloat)(y[i]*y[j]*(this->*kernel_function)(i,j));
  1135. }
  1136. return data;
  1137. }
  1139. void swap_index(int i, int j) const
  1140. {
  1141. cache->swap_index(i,j);
  1142. Kernel::swap_index(i,j);
  1143. swap(y[i],y[j]);
  1144. }
  1146. ~SVC_Q()
  1147. {
  1148. delete[] y;
  1149. delete cache;
  1150. }
  1151. private:
  1152. schar *y;
  1153. Cache *cache;
  1154. };
  1156. class ONE_CLASS_Q: public Kernel
  1157. {
  1158. public:
  1159. ONE_CLASS_Q(const svm_problem& prob, const svm_parameter& param)
  1160. :Kernel(prob.l, prob.x, param)
  1161. {
  1162. cache = new Cache(prob.l,(int)(param.cache_size*(1<<20)));
  1163. }
  1165. Qfloat *get_Q(int i, int len) const
  1166. {
  1167. Qfloat *data;
  1168. int start;
  1169. if((start = cache->get_data(i,&data,len)) < len)
  1170. {
  1171. for(int j=start;j<len;j++)
  1172. data[j] = (Qfloat)(this->*kernel_function)(i,j);
  1173. }
  1174. return data;
  1175. }
  1177. void swap_index(int i, int j) const
  1178. {
  1179. cache->swap_index(i,j);
  1180. Kernel::swap_index(i,j);
  1181. }
  1183. ~ONE_CLASS_Q()
  1184. {
  1185. delete cache;
  1186. }
  1187. private:
  1188. Cache *cache;
  1189. };
  1191. class SVR_Q: public Kernel
  1192. { 
  1193. public:
  1194. SVR_Q(const svm_problem& prob, const svm_parameter& param)
  1195. :Kernel(prob.l, prob.x, param)
  1196. {
  1197. l = prob.l;
  1198. cache = new Cache(l,(int)(param.cache_size*(1<<20)));
  1199. sign = new schar[2*l];
  1200. index = new int[2*l];
  1201. for(int k=0;k<l;k++)
  1202. {
  1203. sign[k] = 1;
  1204. sign[k+l] = -1;
  1205. index[k] = k;
  1206. index[k+l] = k;
  1207. }
  1208. buffer[0] = new Qfloat[2*l];
  1209. buffer[1] = new Qfloat[2*l];
  1210. next_buffer = 0;
  1211. }
  1213. void swap_index(int i, int j) const
  1214. {
  1215. swap(sign[i],sign[j]);
  1216. swap(index[i],index[j]);
  1217. }
  1219. Qfloat *get_Q(int i, int len) const
  1220. {
  1221. Qfloat *data;
  1222. int real_i = index[i];
  1223. if(cache->get_data(real_i,&data,l) < l)
  1224. {
  1225. for(int j=0;j<l;j++)
  1226. data[j] = (Qfloat)(this->*kernel_function)(real_i,j);
  1227. }
  1229. // reorder and copy
  1230. Qfloat *buf = buffer[next_buffer];
  1231. next_buffer = 1 - next_buffer;
  1232. schar si = sign[i];
  1233. for(int j=0;j<len;j++)
  1234. buf[j] = si * sign[j] * data[index[j]];
  1235. return buf;
  1236. }
  1238. ~SVR_Q()
  1239. {
  1240. delete cache;
  1241. delete[] sign;
  1242. delete[] index;
  1243. delete[] buffer[0];
  1244. delete[] buffer[1];
  1245. }
  1246. private:
  1247. int l;
  1248. Cache *cache;
  1249. schar *sign;
  1250. int *index;
  1251. mutable int next_buffer;
  1252. Qfloat* buffer[2];
  1253. };
  1255. //
  1256. // construct and solve various formulations
  1257. //
  1258. static void solve_c_svc(
  1259. const svm_problem *prob, const svm_parameter* param,
  1260. double *alpha, Solver::SolutionInfo* si, double Cp, double Cn)
  1261. {
  1262. int l = prob->l;
  1263. double *minus_ones = new double[l];
  1264. schar *y = new schar[l];
  1266. int i;
  1268. for(i=0;i<l;i++)
  1269. {
  1270. alpha[i] = 0;
  1271. minus_ones[i] = -1;
  1272. if(prob->y[i] > 0) y[i] = +1; else y[i]=-1;
  1273. }
  1275. Solver s;
  1276. s.Solve(l, SVC_Q(*prob,*param,y), minus_ones, y,
  1277. alpha, Cp, Cn, param->eps, si, param->shrinking);
  1279. double sum_alpha=0;
  1280. for(i=0;i<l;i++)
  1281. sum_alpha += alpha[i];
  1283. info("nu = %fn", sum_alpha/(param->C*prob->l));
  1285. for(i=0;i<l;i++)
  1286. alpha[i] *= y[i];
  1288. delete[] minus_ones;
  1289. delete[] y;
  1290. }
  1292. static void solve_nu_svc(
  1293. const svm_problem *prob, const svm_parameter *param,
  1294. double *alpha, Solver::SolutionInfo* si)
  1295. {
  1296. int i;
  1297. int l = prob->l;
  1298. double nu = param->nu;
  1300. int y_pos = 0;
  1301. int y_neg = 0;
  1302. schar *y = new schar[l];
  1304. for(i=0;i<l;i++)
  1305. if(prob->y[i]>0)
  1306. {
  1307. y[i] = +1;
  1308. ++y_pos;
  1309. }
  1310. else
  1311. {
  1312. y[i] = -1;
  1313. ++y_neg;
  1314. }
  1316. if(nu < 0 || nu*l/2 > min(y_pos,y_neg))
  1317. {
  1318. fprintf(stderr,"specified nu is infeasiblen");
  1319. exit(1);
  1320. }
  1322. double sum_pos = nu*l/2;
  1323. double sum_neg = nu*l/2;
  1325. for(i=0;i<l;i++)
  1326. if(y[i] == +1)
  1327. {
  1328. alpha[i] = min(1.0,sum_pos);
  1329. sum_pos -= alpha[i];
  1330. }
  1331. else
  1332. {
  1333. alpha[i] = min(1.0,sum_neg);
  1334. sum_neg -= alpha[i];
  1335. }
  1337. double *zeros = new double[l];
  1339. for(i=0;i<l;i++)
  1340. zeros[i] = 0;
  1342. Solver_NU s;
  1343. s.Solve(l, SVC_Q(*prob,*param,y), zeros, y,
  1344. alpha, 1.0, 1.0, param->eps, si,  param->shrinking);
  1345. double r = si->r;
  1347. info("C = %fn",1/r);
  1349. for(i=0;i<l;i++)
  1350. alpha[i] *= y[i]/r;
  1352. si->rho /= r;
  1353. si->obj /= (r*r);
  1354. si->upper_bound_p = 1/r;
  1355. si->upper_bound_n = 1/r;
  1357. delete[] y;
  1358. delete[] zeros;
  1359. }
  1361. static void solve_one_class(
  1362. const svm_problem *prob, const svm_parameter *param,
  1363. double *alpha, Solver::SolutionInfo* si)
  1364. {
  1365. int l = prob->l;
  1366. double *zeros = new double[l];
  1367. schar *ones = new schar[l];
  1368. int i;
  1370. int n = (int)(param->nu*prob->l); // # of alpha's at upper bound
  1371. if(n>=prob->l)
  1372. {
  1373. fprintf(stderr,"nu must be in (0,1)n");
  1374. exit(1);
  1375. }
  1376. for(i=0;i<n;i++)
  1377. alpha[i] = 1;
  1378. alpha[n] = param->nu * prob->l - n;
  1379. for(i=n+1;i<l;i++)
  1380. alpha[i] = 0;
  1382. for(i=0;i<l;i++)
  1383. {
  1384. zeros[i] = 0;
  1385. ones[i] = 1;
  1386. }
  1388. Solver s;
  1389. s.Solve(l, ONE_CLASS_Q(*prob,*param), zeros, ones,
  1390. alpha, 1.0, 1.0, param->eps, si, param->shrinking);
  1392. delete[] zeros;
  1393. delete[] ones;
  1394. }
  1396. static void solve_epsilon_svr(
  1397. const svm_problem *prob, const svm_parameter *param,
  1398. double *alpha, Solver::SolutionInfo* si)
  1399. {
  1400. int l = prob->l;
  1401. double *alpha2 = new double[2*l];
  1402. double *linear_term = new double[2*l];
  1403. schar *y = new schar[2*l];
  1404. int i;
  1406. for(i=0;i<l;i++)
  1407. {
  1408. alpha2[i] = 0;
  1409. linear_term[i] = param->p - prob->y[i];
  1410. y[i] = 1;
  1412. alpha2[i+l] = 0;
  1413. linear_term[i+l] = param->p + prob->y[i];
  1414. y[i+l] = -1;
  1415. }
  1417. Solver s;
  1418. s.Solve(2*l, SVR_Q(*prob,*param), linear_term, y,
  1419. alpha2, param->C, param->C, param->eps, si, param->shrinking);
  1421. double sum_alpha = 0;
  1422. for(i=0;i<l;i++)
  1423. {
  1424. alpha[i] = alpha2[i] - alpha2[i+l];
  1425. sum_alpha += fabs(alpha[i]);
  1426. }
  1427. info("nu = %fn",sum_alpha/(param->C*l));
  1429. delete[] alpha2;
  1430. delete[] linear_term;
  1431. delete[] y;
  1432. }
  1434. static void solve_nu_svr(
  1435. const svm_problem *prob, const svm_parameter *param,
  1436. double *alpha, Solver::SolutionInfo* si)
  1437. {
  1438. if(param->nu < 0 || param->nu > 1)
  1439. {
  1440. fprintf(stderr,"specified nu is out of rangen");
  1441. exit(1);
  1442. }
  1444. int l = prob->l;
  1445. double C = param->C;
  1446. double *alpha2 = new double[2*l];
  1447. double *linear_term = new double[2*l];
  1448. schar *y = new schar[2*l];
  1449. int i;
  1451. double sum = C * param->nu * l / 2;
  1452. for(i=0;i<l;i++)
  1453. {
  1454. alpha2[i] = alpha2[i+l] = min(sum,C);
  1455. sum -= alpha2[i];
  1457. linear_term[i] = - prob->y[i];
  1458. y[i] = 1;
  1460. linear_term[i+l] = prob->y[i];
  1461. y[i+l] = -1;
  1462. }
  1464. Solver_NU s;
  1465. s.Solve(2*l, SVR_Q(*prob,*param), linear_term, y,
  1466. alpha2, C, C, param->eps, si, param->shrinking);
  1468. info("epsilon = %fn",-si->r);
  1470. for(i=0;i<l;i++)
  1471. alpha[i] = alpha2[i] - alpha2[i+l];
  1473. delete[] alpha2;
  1474. delete[] linear_term;
  1475. delete[] y;
  1476. }
  1478. //
  1479. // decision_function
  1480. //
  1481. struct decision_function
  1482. {
  1483. double *alpha;
  1484. double rho;
  1485. };
  1487. decision_function svm_train_one(
  1488. const svm_problem *prob, const svm_parameter *param,
  1489. double Cp, double Cn)
  1490. {
  1491. double *alpha = Malloc(double,prob->l);
  1492. Solver::SolutionInfo si;
  1493. switch(param->svm_type)
  1494. {
  1495. case C_SVC:
  1496. solve_c_svc(prob,param,alpha,&si,Cp,Cn);
  1497. break;
  1498. case NU_SVC:
  1499. solve_nu_svc(prob,param,alpha,&si);
  1500. break;
  1501. case ONE_CLASS:
  1502. solve_one_class(prob,param,alpha,&si);
  1503. break;
  1504. case EPSILON_SVR:
  1505. solve_epsilon_svr(prob,param,alpha,&si);
  1506. break;
  1507. case NU_SVR:
  1508. solve_nu_svr(prob,param,alpha,&si);
  1509. break;
  1510. }
  1512. info("obj = %f, rho = %fn",si.obj,si.rho);
  1514. // output SVs
  1516. int nSV = 0;
  1517. int nBSV = 0;
  1518. for(int i=0;i<prob->l;i++)
  1519. {
  1520. if(fabs(alpha[i]) > 0)
  1521. {
  1522. ++nSV;
  1523. if(prob->y[i] > 0)
  1524. {
  1525. if(fabs(alpha[i]) >= si.upper_bound_p)
  1526. ++nBSV;
  1527. }
  1528. else
  1529. {
  1530. if(fabs(alpha[i]) >= si.upper_bound_n)
  1531. ++nBSV;
  1532. }
  1533. }
  1534. }
  1536. info("nSV = %d, nBSV = %dn",nSV,nBSV);
  1538. decision_function f;
  1539. f.alpha = alpha;
  1540. f.rho = si.rho;
  1541. return f;
  1542. }
  1544. //
  1545. // svm_model
  1546. //
  1547. struct svm_model
  1548. {
  1549. svm_parameter param; // parameter
  1550. int nr_class; // number of classes, = 2 in regression/one class svm
  1551. int l; // total #SV
  1552. svm_node **SV; // SVs (SV[l])
  1553. double **sv_coef; // coefficients for SVs in decision functions (sv_coef[n-1][l])
  1554. double *rho; // constants in decision functions (rho[n*(n-1)/2])
  1556. // for classification only
  1558. int *label; // label of each class (label[n])
  1559. int *nSV; // number of SVs for each class (nSV[n])
  1560. // nSV[0] + nSV[1] + ... + nSV[n-1] = l
  1561. // XXX
  1562. int free_sv; // 1 if svm_model is created by svm_load_model
  1563. // 0 if svm_model is created by svm_train
  1564. };
  1566. //
  1567. // Interface functions
  1568. //
  1569. svm_model *svm_train(const svm_problem *prob, const svm_parameter *param)
  1570. {
  1571. svm_model *model = Malloc(svm_model,1);
  1572. model->param = *param;
  1573. model->free_sv = 0; // XXX
  1575. if(param->svm_type == ONE_CLASS ||
  1576.    param->svm_type == EPSILON_SVR ||
  1577.    param->svm_type == NU_SVR)
  1578. {
  1579. // regression or one-class-svm
  1580. model->nr_class = 2;
  1581. model->label = NULL;
  1582. model->nSV = NULL;
  1583. model->sv_coef = Malloc(double *,1);
  1584. decision_function f = svm_train_one(prob,param,0,0);
  1585. model->rho = Malloc(double,1);
  1586. model->rho[0] = f.rho;
  1588. int nSV = 0;
  1589. int i;
  1590. for(i=0;i<prob->l;i++)
  1591. if(fabs(f.alpha[i]) > 0) ++nSV;
  1592. model->l = nSV;
  1593. model->SV = Malloc(svm_node *,nSV);
  1594. model->sv_coef[0] = Malloc(double,nSV);
  1595. int j = 0;
  1596. for(i=0;i<prob->l;i++)
  1597. if(fabs(f.alpha[i]) > 0)
  1598. {
  1599. model->SV[j] = prob->x[i];
  1600. model->sv_coef[0][j] = f.alpha[i];
  1601. ++j;
  1602. }
  1604. free(f.alpha);
  1605. }
  1606. else
  1607. {
  1608. // classification
  1609. // find out the number of classes
  1610. int l = prob->l;
  1611. int max_nr_class = 16;
  1612. int nr_class = 0;
  1613. int *label = Malloc(int,max_nr_class);
  1614. int *count = Malloc(int,max_nr_class);
  1615. int *index = Malloc(int,l);
  1617. int i;
  1618. for(i=0;i<l;i++)
  1619. {
  1620. int this_label = (int)prob->y[i];
  1621. int j;
  1622. for(j=0;j<nr_class;j++)
  1623. if(this_label == label[j])
  1624. {
  1625. ++count[j];
  1626. break;
  1627. }
  1628. index[i] = j;
  1629. if(j == nr_class)
  1630. {
  1631. if(nr_class == max_nr_class)
  1632. {
  1633. max_nr_class *= 2;
  1634. label = (int *)realloc(label,max_nr_class*sizeof(int));
  1635. count = (int *)realloc(count,max_nr_class*sizeof(int));
  1636. }
  1637. label[nr_class] = this_label;
  1638. count[nr_class] = 1;
  1639. ++nr_class;
  1640. }
  1641. }
  1643. // group training data of the same class
  1645. int *start = Malloc(int,nr_class);
  1646. start[0] = 0;
  1647. for(i=1;i<nr_class;i++)
  1648. start[i] = start[i-1]+count[i-1];
  1650. svm_node **x = Malloc(svm_node *,l);
  1652. for(i=0;i<l;i++)
  1653. {
  1654. x[start[index[i]]] = prob->x[i];
  1655. ++start[index[i]];
  1656. }
  1658. start[0] = 0;
  1659. for(i=1;i<nr_class;i++)
  1660. start[i] = start[i-1]+count[i-1];
  1662. // calculate weighted C
  1664. double *weighted_C = Malloc(double, nr_class);
  1665. for(i=0;i<nr_class;i++)
  1666. weighted_C[i] = param->C;
  1667. for(i=0;i<param->nr_weight;i++)
  1668. {
  1669. int j;
  1670. for(j=0;j<nr_class;j++)
  1671. if(param->weight_label[i] == label[j])
  1672. break;
  1673. if(j == nr_class)
  1674. fprintf(stderr,"warning: class label %d specified in weight is not foundn", param->weight_label[i]);
  1675. else
  1676. weighted_C[j] *= param->weight[i];
  1677. }
  1679. // train n*(n-1)/2 models
  1681. bool *nonzero = Malloc(bool,l);
  1682. for(i=0;i<l;i++)
  1683. nonzero[i] = false;
  1684. decision_function *f = Malloc(decision_function,nr_class*(nr_class-1)/2);
  1686. int p = 0;
  1687. for(i=0;i<nr_class;i++)
  1688. for(int j=i+1;j<nr_class;j++)
  1689. {
  1690. svm_problem sub_prob;
  1691. int si = start[i], sj = start[j];
  1692. int ci = count[i], cj = count[j];
  1693. sub_prob.l = ci+cj;
  1694. sub_prob.x = Malloc(svm_node *,sub_prob.l);
  1695. sub_prob.y = Malloc(double,sub_prob.l);
  1696. int k;
  1697. for(k=0;k<ci;k++)
  1698. {
  1699. sub_prob.x[k] = x[si+k];
  1700. sub_prob.y[k] = +1;
  1701. }
  1702. for(k=0;k<cj;k++)
  1703. {
  1704. sub_prob.x[ci+k] = x[sj+k];
  1705. sub_prob.y[ci+k] = -1;
  1706. }
  1708. f[p] = svm_train_one(&sub_prob,param,weighted_C[i],weighted_C[j]);
  1709. for(k=0;k<ci;k++)
  1710. if(!nonzero[si+k] && fabs(f[p].alpha[k]) > 0)
  1711. nonzero[si+k] = true;
  1712. for(k=0;k<cj;k++)
  1713. if(!nonzero[sj+k] && fabs(f[p].alpha[ci+k]) > 0)
  1714. nonzero[sj+k] = true;
  1715. free(sub_prob.x);
  1716. free(sub_prob.y);
  1717. ++p;
  1718. }
  1720. // build output
  1722. model->nr_class = nr_class;
  1724. model->label = Malloc(int,nr_class);
  1725. for(i=0;i<nr_class;i++)
  1726. model->label[i] = label[i];
  1728. model->rho = Malloc(double,nr_class*(nr_class-1)/2);
  1729. for(i=0;i<nr_class*(nr_class-1)/2;i++)
  1730. model->rho[i] = f[i].rho;
  1732. int total_sv = 0;
  1733. int *nz_count = Malloc(int,nr_class);
  1734. model->nSV = Malloc(int,nr_class);
  1735. for(i=0;i<nr_class;i++)
  1736. {
  1737. int nSV = 0;
  1738. for(int j=0;j<count[i];j++)
  1739. if(nonzero[start[i]+j])
  1740. {
  1741. ++nSV;
  1742. ++total_sv;
  1743. }
  1744. model->nSV[i] = nSV;
  1745. nz_count[i] = nSV;
  1746. }
  1748. info("Total nSV = %dn",total_sv);
  1750. model->l = total_sv;
  1751. model->SV = Malloc(svm_node *,total_sv);
  1752. p = 0;
  1753. for(i=0;i<l;i++)
  1754. if(nonzero[i]) model->SV[p++] = x[i];
  1756. int *nz_start = Malloc(int,nr_class);
  1757. nz_start[0] = 0;
  1758. for(i=1;i<nr_class;i++)
  1759. nz_start[i] = nz_start[i-1]+nz_count[i-1];
  1761. model->sv_coef = Malloc(double *,nr_class-1);
  1762. for(i=0;i<nr_class-1;i++)
  1763. model->sv_coef[i] = Malloc(double,total_sv);
  1765. p = 0;
  1766. for(i=0;i<nr_class;i++)
  1767. for(int j=i+1;j<nr_class;j++)
  1768. {
  1769. // classifier (i,j): coefficients with
  1770. // i are in sv_coef[j-1][nz_start[i]...],
  1771. // j are in sv_coef[i][nz_start[j]...]
  1773. int si = start[i];
  1774. int sj = start[j];
  1775. int ci = count[i];
  1776. int cj = count[j];
  1778. int q = nz_start[i];
  1779. int k;
  1780. for(k=0;k<ci;k++)
  1781. if(nonzero[si+k])
  1782. model->sv_coef[j-1][q++] = f[p].alpha[k];
  1783. q = nz_start[j];
  1784. for(k=0;k<cj;k++)
  1785. if(nonzero[sj+k])
  1786. model->sv_coef[i][q++] = f[p].alpha[ci+k];
  1787. ++p;
  1788. }
  1790. free(label);
  1791. free(count);
  1792. free(index);
  1793. free(start);
  1794. free(x);
  1795. free(weighted_C);
  1796. free(nonzero);
  1797. for(i=0;i<nr_class*(nr_class-1)/2;i++)
  1798. free(f[i].alpha);
  1799. free(f);
  1800. free(nz_count);
  1801. free(nz_start);
  1802. }
  1803. return model;
  1804. }
  1806. double svm_predict(const svm_model *model, const svm_node *x)
  1807. {
  1808. if(model->param.svm_type == ONE_CLASS ||
  1809.    model->param.svm_type == EPSILON_SVR ||
  1810.    model->param.svm_type == NU_SVR)
  1811. {
  1812. double *sv_coef = model->sv_coef[0];
  1813. double sum = 0;
  1814. for(int i=0;i<model->l;i++)
  1815. sum += sv_coef[i] * Kernel::k_function(x,model->SV[i],model->param);
  1816. sum -= model->rho[0];
  1817. if(model->param.svm_type == ONE_CLASS)
  1818. return (sum>0)?1:-1;
  1819. else
  1820. return sum;
  1821. }
  1822. else
  1823. {
  1824. int i;
  1825. int nr_class = model->nr_class;
  1826. int l = model->l;
  1828. double *kvalue = Malloc(double,l);
  1829. for(i=0;i<l;i++)
  1830. kvalue[i] = Kernel::k_function(x,model->SV[i],model->param);
  1832. int *start = Malloc(int,nr_class);
  1833. start[0] = 0;
  1834. for(i=1;i<nr_class;i++)
  1835. start[i] = start[i-1]+model->nSV[i-1];
  1837. int *vote = Malloc(int,nr_class);
  1838. for(i=0;i<nr_class;i++)
  1839. vote[i] = 0;
  1840. int p=0;
  1841. for(i=0;i<nr_class;i++)
  1842. for(int j=i+1;j<nr_class;j++)
  1843. {
  1844. double sum = 0;
  1845. int si = start[i];
  1846. int sj = start[j];
  1847. int ci = model->nSV[i];
  1848. int cj = model->nSV[j];
  1850. int k;
  1851. double *coef1 = model->sv_coef[j-1];
  1852. double *coef2 = model->sv_coef[i];
  1853. for(k=0;k<ci;k++)
  1854. sum += coef1[si+k] * kvalue[si+k];
  1855. for(k=0;k<cj;k++)
  1856. sum += coef2[sj+k] * kvalue[sj+k];
  1857. sum -= model->rho[p++];
  1858. if(sum > 0)
  1859. ++vote[i];
  1860. else
  1861. ++vote[j];
  1862. }
  1864. int vote_max_idx = 0;
  1865. for(i=1;i<nr_class;i++)
  1866. if(vote[i] > vote[vote_max_idx])
  1867. vote_max_idx = i;
  1868. free(kvalue);
  1869. free(start);
  1870. free(vote);
  1871. return model->label[vote_max_idx];
  1872. }
  1873. }
  1875. const char *svm_type_table[] =
  1876. {
  1877. "c_svc","nu_svc","one_class","epsilon_svr","nu_svr",NULL
  1878. };
  1880. const char *kernel_type_table[]=
  1881. {
  1882. "linear","polynomial","rbf","sigmoid",NULL
  1883. };
  1885. int svm_save_model(const char *model_file_name, const svm_model *model)
  1886. {
  1887. FILE *fp = fopen(model_file_name,"w");
  1888. if(fp==NULL) return -1;
  1890. const svm_parameter& param = model->param;
  1892. fprintf(fp,"svm_type %sn", svm_type_table[param.svm_type]);
  1893. fprintf(fp,"kernel_type %sn", kernel_type_table[param.kernel_type]);
  1895. if(param.kernel_type == POLY)
  1896. fprintf(fp,"degree %gn", param.degree);
  1898. if(param.kernel_type == POLY || param.kernel_type == RBF || param.kernel_type == SIGMOID)
  1899. fprintf(fp,"gamma %gn", param.gamma);
  1901. if(param.kernel_type == POLY || param.kernel_type == SIGMOID)
  1902. fprintf(fp,"coef0 %gn", param.coef0);
  1904. int nr_class = model->nr_class;
  1905. int l = model->l;
  1906. fprintf(fp, "nr_class %dn", nr_class);
  1907. fprintf(fp, "total_sv %dn",l);
  1909. {
  1910. fprintf(fp, "rho");
  1911. for(int i=0;i<nr_class*(nr_class-1)/2;i++)
  1912. fprintf(fp," %g",model->rho[i]);
  1913. fprintf(fp, "n");
  1914. }
  1916. if(model->label)
  1917. {
  1918. fprintf(fp, "label");
  1919. for(int i=0;i<nr_class;i++)
  1920. fprintf(fp," %d",model->label[i]);
  1921. fprintf(fp, "n");
  1922. }
  1924. if(model->nSV)
  1925. {
  1926. fprintf(fp, "nr_sv");
  1927. for(int i=0;i<nr_class;i++)
  1928. fprintf(fp," %d",model->nSV[i]);
  1929. fprintf(fp, "n");
  1930. }
  1932. fprintf(fp, "SVn");
  1933. const double * const *sv_coef = model->sv_coef;
  1934. const svm_node * const *SV = model->SV;
  1936. for(int i=0;i<l;i++)
  1937. {
  1938. for(int j=0;j<nr_class-1;j++)
  1939. fprintf(fp, "%.16g ",sv_coef[j][i]);
  1941. const svm_node *p = SV[i];
  1942. while(p->index != -1)
  1943. {
  1944. fprintf(fp,"%d:%.8g ",p->index,p->value);
  1945. p++;
  1946. }
  1947. fprintf(fp, "n");
  1948. }
  1950. fclose(fp);
  1951. return 0;
  1952. }
  1954. svm_model *svm_load_model(const char *model_file_name)
  1955. {
  1956. FILE *fp = fopen(model_file_name,"rb");
  1957. if(fp==NULL) return NULL;
  1959. // read parameters
  1961. svm_model *model = (svm_model *)malloc(sizeof(svm_model));
  1962. svm_parameter& param = model->param;
  1963. model->label = NULL;
  1964. model->nSV = NULL;
  1966. char cmd[81];
  1967. while(1)
  1968. {
  1969. fscanf(fp,"%80s",cmd);
  1971. if(strcmp(cmd,"svm_type")==0)
  1972. {
  1973. fscanf(fp,"%80s",cmd);
  1974. int i;
  1975. for(i=0;svm_type_table[i];i++)
  1976. {
  1977. if(strcmp(svm_type_table[i],cmd)==0)
  1978. {
  1979. param.svm_type=i;
  1980. break;
  1981. }
  1982. }
  1983. if(svm_type_table[i] == NULL)
  1984. {
  1985. fprintf(stderr,"unknown svm type.n");
  1986. exit(1);
  1987. }
  1988. }
  1989. else if(strcmp(cmd,"kernel_type")==0)
  1990. {
  1991. fscanf(fp,"%80s",cmd);
  1992. int i;
  1993. for(i=0;kernel_type_table[i];i++)
  1994. {
  1995. if(strcmp(kernel_type_table[i],cmd)==0)
  1996. {
  1997. param.kernel_type=i;
  1998. break;
  1999. }
  2000. }
  2001. if(kernel_type_table[i] == NULL)
  2002. {
  2003. fprintf(stderr,"unknown kernel function.n");
  2004. exit(1);
  2005. }
  2006. }
  2007. else if(strcmp(cmd,"degree")==0)
  2008. fscanf(fp,"%lf",&param.degree);
  2009. else if(strcmp(cmd,"gamma")==0)
  2010. fscanf(fp,"%lf",&param.gamma);
  2011. else if(strcmp(cmd,"coef0")==0)
  2012. fscanf(fp,"%lf",&param.coef0);
  2013. else if(strcmp(cmd,"nr_class")==0)
  2014. fscanf(fp,"%d",&model->nr_class);
  2015. else if(strcmp(cmd,"total_sv")==0)
  2016. fscanf(fp,"%d",&model->l);
  2017. else if(strcmp(cmd,"rho")==0)
  2018. {
  2019. int n = model->nr_class * (model->nr_class-1)/2;
  2020. model->rho = Malloc(double,n);
  2021. for(int i=0;i<n;i++)
  2022. fscanf(fp,"%lf",&model->rho[i]);
  2023. }
  2024. else if(strcmp(cmd,"label")==0)
  2025. {
  2026. int n = model->nr_class;
  2027. model->label = Malloc(int,n);
  2028. for(int i=0;i<n;i++)
  2029. fscanf(fp,"%d",&model->label[i]);
  2030. }
  2031. else if(strcmp(cmd,"nr_sv")==0)
  2032. {
  2033. int n = model->nr_class;
  2034. model->nSV = Malloc(int,n);
  2035. for(int i=0;i<n;i++)
  2036. fscanf(fp,"%d",&model->nSV[i]);
  2037. }
  2038. else if(strcmp(cmd,"SV")==0)
  2039. {
  2040. while(1)
  2041. {
  2042. int c = getc(fp);
  2043. if(c==EOF || c=='n') break;
  2044. }
  2045. break;
  2046. }
  2047. else
  2048. {
  2049. fprintf(stderr,"unknown text in model filen");
  2050. exit(1);
  2051. }
  2052. }
  2054. // read sv_coef and SV
  2056. int elements = 0;
  2057. long pos = ftell(fp);
  2059. while(1)
  2060. {
  2061. int c = fgetc(fp);
  2062. switch(c)
  2063. {
  2064. case 'n':
  2065. // count the '-1' element
  2066. case ':':
  2067. ++elements;
  2068. break;
  2069. case EOF:
  2070. goto out;
  2071. default:
  2072. ;
  2073. }
  2074. }
  2075. out:
  2076. fseek(fp,pos,SEEK_SET);
  2078. int m = model->nr_class - 1;
  2079. int l = model->l;
  2080. model->sv_coef = Malloc(double *,m);
  2081. int i;
  2082. for(i=0;i<m;i++)
  2083. model->sv_coef[i] = Malloc(double,l);
  2084. model->SV = Malloc(svm_node*,l);
  2085. svm_node *x_space = Malloc(svm_node,elements);
  2087. int j=0;
  2088. for(i=0;i<l;i++)
  2089. {
  2090. model->SV[i] = &x_space[j];
  2091. for(int k=0;k<m;k++)
  2092. fscanf(fp,"%lf",&model->sv_coef[k][i]);
  2093. while(1)
  2094. {
  2095. int c;
  2096. do {
  2097. c = getc(fp);
  2098. if(c=='n') goto out2;
  2099. } while(isspace(c));
  2100. ungetc(c,fp);
  2101. fscanf(fp,"%d:%lf",&(x_space[j].index),&(x_space[j].value));
  2102. ++j;
  2103. }
  2104. out2:
  2105. x_space[j++].index = -1;
  2106. }
  2108. fclose(fp);
  2110. model->free_sv = 1; // XXX
  2111. return model;
  2112. }
  2114. void svm_destroy_model(svm_model* model)
  2115. {
  2116. if(model->free_sv)
  2117. free((void *)(model->SV[0]));
  2118. for(int i=0;i<model->nr_class-1;i++)
  2119. free(model->sv_coef[i]);
  2120. free(model->SV);
  2121. free(model->sv_coef);
  2122. free(model->rho);
  2123. free(model->label);
  2124. free(model->nSV);
  2125. free(model);
  2126. }