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

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

.net编程

开发平台：

Java

svm.m4：源码内容
  1. define(`swap',`do {$1 _=$2; $2=$3; $3=_;} while(false)')
  2. define(`Qfloat',`float')
  3. define(`SIZE_OF_QFLOAT',4)
  4. package libsvm;
  5. import java.io.*;
  6. import java.util.*;
  8. //
  9. // Kernel Cache
  10. //
  11. // l is the number of total data items
  12. // size is the cache size limit in bytes
  13. //
  14. class Cache {
  15. private final int l;
  16. private int size;
  17. private final class head_t
  18. {
  19. head_t prev, next; // a cicular list
  20. Qfloat[] data;
  21. int len; // data[0,len) is cached in this entry
  22. }
  23. private final head_t[] head;
  24. private head_t lru_head;
  26. Cache(int l_, int size_)
  27. {
  28. l = l_;
  29. size = size_;
  30. head = new head_t[l];
  31. for(int i=0;i<l;i++) head[i] = new head_t();
  32. size /= SIZE_OF_QFLOAT;
  33. size -= l * (16/SIZE_OF_QFLOAT); // sizeof(head_t) == 16
  34. lru_head = new head_t();
  35. lru_head.next = lru_head.prev = lru_head;
  36. }
  38. private void lru_delete(head_t h)
  39. {
  40. // delete from current location
  41. h.prev.next = h.next;
  42. h.next.prev = h.prev;
  43. }
  45. private void lru_insert(head_t h)
  46. {
  47. // insert to last position
  48. h.next = lru_head;
  49. h.prev = lru_head.prev;
  50. h.prev.next = h;
  51. h.next.prev = h;
  52. }
  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. // java: simulate pointer using single-element array
  58. int get_data(int index, Qfloat[][] data, int len)
  59. {
  60. head_t h = head[index];
  61. if(h.len > 0) lru_delete(h);
  62. int more = len - h.len;
  64. if(more > 0)
  65. {
  66. // free old space
  67. while(size < more)
  68. {
  69. head_t old = lru_head.next;
  70. lru_delete(old);
  71. size += old.len;
  72. old.data = null;
  73. old.len = 0;
  74. }
  76. // allocate new space
  77. Qfloat[] new_data = new Qfloat[len];
  78. if(h.data != null) System.arraycopy(h.data,0,new_data,0,h.len);
  79. h.data = new_data;
  80. size -= more;
  81. swap(int,h.len,len);
  82. }
  84. lru_insert(h);
  85. data[0] = h.data;
  86. return len;
  87. }
  89. void swap_index(int i, int j)
  90. {
  91. if(i==j) return;
  93. if(head[i].len > 0) lru_delete(head[i]);
  94. if(head[j].len > 0) lru_delete(head[j]);
  95. swap(Qfloat[],head[i].data,head[j].data);
  96. swap(int,head[i].len,head[j].len);
  97. if(head[i].len > 0) lru_insert(head[i]);
  98. if(head[j].len > 0) lru_insert(head[j]);
  100. if(i>j) swap(int,i,j);
  101. for(head_t h = lru_head.next; h!=lru_head; h=h.next)
  102. {
  103. if(h.len > i)
  104. {
  105. if(h.len > j)
  106. swap(Qfloat,h.data[i],h.data[j]);
  107. else
  108. {
  109. // give up
  110. lru_delete(h);
  111. size += h.len;
  112. h.data = null;
  113. h.len = 0;
  114. }
  115. }
  116. }
  117. }
  118. }
  120. //
  121. // Kernel evaluation
  122. //
  123. // the static method k_function is for doing single kernel evaluation
  124. // the constructor of Kernel prepares to calculate the l*l kernel matrix
  125. // the member function get_Q is for getting one column from the Q Matrix
  126. //
  127. abstract class Kernel {
  128. private svm_node[][] x;
  129. private final double[] x_square;
  131. // svm_parameter
  132. private final int kernel_type;
  133. private final double degree;
  134. private final double gamma;
  135. private final double coef0;
  137. abstract Qfloat[] get_Q(int column, int len);
  139. void swap_index(int i, int j)
  140. {
  141. swap(svm_node[],x[i],x[j]);
  142. if(x_square != null) swap(double,x_square[i],x_square[j]);
  143. }
  145. private static double tanh(double x)
  146. {
  147. double e = Math.exp(x);
  148. return 1.0-2.0/(e*e+1);
  149. }
  151. double kernel_function(int i, int j)
  152. {
  153. switch(kernel_type)
  154. {
  155. case svm_parameter.LINEAR:
  156. return dot(x[i],x[j]);
  157. case svm_parameter.POLY:
  158. return Math.pow(gamma*dot(x[i],x[j])+coef0,degree);
  159. case svm_parameter.RBF:
  160. return Math.exp(-gamma*(x_square[i]+x_square[j]-2*dot(x[i],x[j])));
  161. case svm_parameter.SIGMOID:
  162. return tanh(gamma*dot(x[i],x[j])+coef0);
  163. default:
  164. System.err.print("unknown kernel function.n");
  165. System.exit(1);
  166. return 0; // java
  167. }
  168. }
  170. Kernel(int l, svm_node[][] x_, svm_parameter param)
  171. {
  172. this.kernel_type = param.kernel_type;
  173. this.degree = param.degree;
  174. this.gamma = param.gamma;
  175. this.coef0 = param.coef0;
  177. x = (svm_node[][])x_.clone();
  179. if(kernel_type == svm_parameter.RBF)
  180. {
  181. x_square = new double[l];
  182. for(int i=0;i<l;i++)
  183. x_square[i] = dot(x[i],x[i]);
  184. }
  185. else x_square = null;
  186. }
  188. static double dot(svm_node[] x, svm_node[] y)
  189. {
  190. double sum = 0;
  191. int xlen = x.length;
  192. int ylen = y.length;
  193. int i = 0;
  194. int j = 0;
  195. while(i < xlen && j < ylen)
  196. {
  197. if(x[i].index == y[j].index)
  198. sum += x[i++].value * y[j++].value;
  199. else
  200. {
  201. if(x[i].index > y[j].index)
  202. ++j;
  203. else
  204. ++i;
  205. }
  206. }
  207. return sum;
  208. }
  210. static double k_function(svm_node[] x, svm_node[] y,
  211. svm_parameter param)
  212. {
  213. switch(param.kernel_type)
  214. {
  215. case svm_parameter.LINEAR:
  216. return dot(x,y);
  217. case svm_parameter.POLY:
  218. return Math.pow(param.gamma*dot(x,y)+param.coef0,param.degree);
  219. case svm_parameter.RBF:
  220. {
  221. double sum = 0;
  222. int xlen = x.length;
  223. int ylen = y.length;
  224. int i = 0;
  225. int j = 0;
  226. while(i < xlen && j < ylen)
  227. {
  228. if(x[i].index == y[j].index)
  229. {
  230. double d = x[i++].value - y[j++].value;
  231. sum += d*d;
  232. }
  233. else if(x[i].index > y[j].index)
  234. {
  235. sum += y[j].value * y[j].value;
  236. ++j;
  237. }
  238. else
  239. {
  240. sum += x[i].value * x[i].value;
  241. ++i;
  242. }
  243. }
  245. while(i < xlen)
  246. {
  247. sum += x[i].value * x[i].value;
  248. ++i;
  249. }
  251. while(j < ylen)
  252. {
  253. sum += y[j].value * y[j].value;
  254. ++j;
  255. }
  257. return Math.exp(-param.gamma*sum);
  258. }
  259. case svm_parameter.SIGMOID:
  260. return tanh(param.gamma*dot(x,y)+param.coef0);
  261. default:
  262. System.err.print("unknown kernel function.n");
  263. System.exit(1);
  264. return 0; // java
  265. }
  266. }
  267. }
  269. // Generalized SMO+SVMlight algorithm
  270. // Solves:
  271. //
  272. // min 0.5(alpha^T Q alpha) + b^T alpha
  273. //
  274. // y^T alpha = delta
  275. // y_i = +1 or -1
  276. // 0 <= alpha_i <= Cp for y_i = 1
  277. // 0 <= alpha_i <= Cn for y_i = -1
  278. //
  279. // Given:
  280. //
  281. // Q, b, y, Cp, Cn, and an initial feasible point alpha
  282. // l is the size of vectors and matrices
  283. // eps is the stopping criterion
  284. //
  285. // solution will be put in alpha, objective value will be put in obj
  286. //
  287. class Solver {
  288. int active_size;
  289. byte[] y;
  290. double[] G; // gradient of objective function
  291. static final byte LOWER_BOUND = 0;
  292. static final byte UPPER_BOUND = 1;
  293. static final byte FREE = 2;
  294. byte[] alpha_status; // LOWER_BOUND, UPPER_BOUND, FREE
  295. double[] alpha;
  296. Kernel Q;
  297. double eps;
  298. double Cp,Cn;
  299. double[] b;
  300. int[] active_set;
  301. double[] G_bar; // gradient, if we treat free variables as 0
  302. int l;
  303. boolean unshrinked; // XXX
  305. static final double INF = java.lang.Double.POSITIVE_INFINITY;
  307. double get_C(int i)
  308. {
  309. return (y[i] > 0)? Cp : Cn;
  310. }
  311. void update_alpha_status(int i)
  312. {
  313. if(alpha[i] >= get_C(i))
  314. alpha_status[i] = UPPER_BOUND;
  315. else if(alpha[i] <= 0)
  316. alpha_status[i] = LOWER_BOUND;
  317. else alpha_status[i] = FREE;
  318. }
  319. boolean is_upper_bound(int i) { return alpha_status[i] == UPPER_BOUND; }
  320. boolean is_lower_bound(int i) { return alpha_status[i] == LOWER_BOUND; }
  321. boolean is_free(int i) {  return alpha_status[i] == FREE; }
  323. // java: information about solution except alpha,
  324. // because we cannot return multiple values otherwise...
  325. static class SolutionInfo {
  326. double obj;
  327. double rho;
  328. double upper_bound_p;
  329. double upper_bound_n;
  330. double r; // for Solver_NU
  331. }
  333. void swap_index(int i, int j)
  334. {
  335. Q.swap_index(i,j);
  336. swap(byte, y[i],y[j]);
  337. swap(double, G[i],G[j]);
  338. swap(byte, alpha_status[i],alpha_status[j]);
  339. swap(double, alpha[i],alpha[j]);
  340. swap(double, b[i],b[j]);
  341. swap(int, active_set[i],active_set[j]);
  342. swap(double, G_bar[i],G_bar[j]);
  343. }
  345. void reconstruct_gradient()
  346. {
  347. // reconstruct inactive elements of G from G_bar and free variables
  349. if(active_size == l) return;
  351. int i;
  352. for(i=active_size;i<l;i++)
  353. G[i] = G_bar[i] + b[i];
  355. for(i=0;i<active_size;i++)
  356. if(is_free(i))
  357. {
  358. Qfloat[] Q_i = Q.get_Q(i,l);
  359. double alpha_i = alpha[i];
  360. for(int j=active_size;j<l;j++)
  361. G[j] += alpha_i * Q_i[j];
  362. }
  363. }
  365. void Solve(int l, Kernel Q, double[] b_, byte[] y_,
  366.    double[] alpha_, double Cp, double Cn, double eps, SolutionInfo si, int shrinking)
  367. {
  368. this.l = l;
  369. this.Q = Q;
  370. b = (double[])b_.clone();
  371. y = (byte[])y_.clone();
  372. alpha = (double[])alpha_.clone();
  373. this.Cp = Cp;
  374. this.Cn = Cn;
  375. this.eps = eps;
  376. this.unshrinked = false;
  378. // initialize alpha_status
  379. {
  380. alpha_status = new byte[l];
  381. for(int i=0;i<l;i++)
  382. update_alpha_status(i);
  383. }
  385. // initialize active set (for shrinking)
  386. {
  387. active_set = new int[l];
  388. for(int i=0;i<l;i++)
  389. active_set[i] = i;
  390. active_size = l;
  391. }
  393. // initialize gradient
  394. {
  395. G = new double[l];
  396. G_bar = new double[l];
  397. int i;
  398. for(i=0;i<l;i++)
  399. {
  400. G[i] = b[i];
  401. G_bar[i] = 0;
  402. }
  403. for(i=0;i<l;i++)
  404. if(!is_lower_bound(i))
  405. {
  406. Qfloat[] Q_i = Q.get_Q(i,l);
  407. double alpha_i = alpha[i];
  408. int j;
  409. for(j=0;j<l;j++)
  410. G[j] += alpha_i*Q_i[j];
  411. if(is_upper_bound(i))
  412. for(j=0;j<l;j++)
  413. G_bar[j] += get_C(i) * Q_i[j];
  414. }
  415. }
  417. // optimization step
  419. int iter = 0;
  420. int counter = Math.min(l,1000)+1;
  421. int[] working_set = new int[2];
  423. while(true)
  424. {
  425. // show progress and do shrinking
  427. if(--counter == 0)
  428. {
  429. counter = Math.min(l,1000);
  430. if(shrinking!=0) do_shrinking();
  431. System.err.print(".");
  432. }
  434. if(select_working_set(working_set)!=0)
  435. {
  436. // reconstruct the whole gradient
  437. reconstruct_gradient();
  438. // reset active set size and check
  439. active_size = l;
  440. System.err.print("*");
  441. if(select_working_set(working_set)!=0)
  442. break;
  443. else
  444. counter = 1; // do shrinking next iteration
  445. }
  447. int i = working_set[0];
  448. int j = working_set[1];
  450. ++iter;
  452. // update alpha[i] and alpha[j], handle bounds carefully
  454. Qfloat[] Q_i = Q.get_Q(i,active_size);
  455. Qfloat[] Q_j = Q.get_Q(j,active_size);
  457. double C_i = get_C(i);
  458. double C_j = get_C(j);
  460. double old_alpha_i = alpha[i];
  461. double old_alpha_j = alpha[j];
  463. if(y[i]!=y[j])
  464. {
  465. double delta = (-G[i]-G[j])/(Q_i[i]+Q_j[j]+2*Q_i[j]);
  466. double diff = alpha[i] - alpha[j];
  467. alpha[i] += delta;
  468. alpha[j] += delta;
  470. if(diff > 0)
  471. {
  472. if(alpha[j] < 0)
  473. {
  474. alpha[j] = 0;
  475. alpha[i] = diff;
  476. }
  477. }
  478. else
  479. {
  480. if(alpha[i] < 0)
  481. {
  482. alpha[i] = 0;
  483. alpha[j] = -diff;
  484. }
  485. }
  486. if(diff > C_i - C_j)
  487. {
  488. if(alpha[i] > C_i)
  489. {
  490. alpha[i] = C_i;
  491. alpha[j] = C_i - diff;
  492. }
  493. }
  494. else
  495. {
  496. if(alpha[j] > C_j)
  497. {
  498. alpha[j] = C_j;
  499. alpha[i] = C_j + diff;
  500. }
  501. }
  502. }
  503. else
  504. {
  505. double delta = (G[i]-G[j])/(Q_i[i]+Q_j[j]-2*Q_i[j]);
  506. double sum = alpha[i] + alpha[j];
  507. alpha[i] -= delta;
  508. alpha[j] += delta;
  509. if(sum > C_i)
  510. {
  511. if(alpha[i] > C_i)
  512. {
  513. alpha[i] = C_i;
  514. alpha[j] = sum - C_i;
  515. }
  516. }
  517. else
  518. {
  519. if(alpha[j] < 0)
  520. {
  521. alpha[j] = 0;
  522. alpha[i] = sum;
  523. }
  524. }
  525. if(sum > C_j)
  526. {
  527. if(alpha[j] > C_j)
  528. {
  529. alpha[j] = C_j;
  530. alpha[i] = sum - C_j;
  531. }
  532. }
  533. else
  534. {
  535. if(alpha[i] < 0)
  536. {
  537. alpha[i] = 0;
  538. alpha[j] = sum;
  539. }
  540. }
  541. }
  543. // update G
  545. double delta_alpha_i = alpha[i] - old_alpha_i;
  546. double delta_alpha_j = alpha[j] - old_alpha_j;
  548. for(int k=0;k<active_size;k++)
  549. {
  550. G[k] += Q_i[k]*delta_alpha_i + Q_j[k]*delta_alpha_j;
  551. }
  553. // update alpha_status and G_bar
  555. {
  556. boolean ui = is_upper_bound(i);
  557. boolean uj = is_upper_bound(j);
  558. update_alpha_status(i);
  559. update_alpha_status(j);
  560. int k;
  561. if(ui != is_upper_bound(i))
  562. {
  563. Q_i = Q.get_Q(i,l);
  564. if(ui)
  565. for(k=0;k<l;k++)
  566. G_bar[k] -= C_i * Q_i[k];
  567. else
  568. for(k=0;k<l;k++)
  569. G_bar[k] += C_i * Q_i[k];
  570. }
  572. if(uj != is_upper_bound(j))
  573. {
  574. Q_j = Q.get_Q(j,l);
  575. if(uj)
  576. for(k=0;k<l;k++)
  577. G_bar[k] -= C_j * Q_j[k];
  578. else
  579. for(k=0;k<l;k++)
  580. G_bar[k] += C_j * Q_j[k];
  581. }
  582. }
  584. }
  586. // calculate rho
  588. si.rho = calculate_rho();
  590. // calculate objective value
  591. {
  592. double v = 0;
  593. int i;
  594. for(i=0;i<l;i++)
  595. v += alpha[i] * (G[i] + b[i]);
  597. si.obj = v/2;
  598. }
  600. // put back the solution
  601. {
  602. for(int i=0;i<l;i++)
  603. alpha_[active_set[i]] = alpha[i];
  604. }
  606. si.upper_bound_p = Cp;
  607. si.upper_bound_n = Cn;
  609. System.out.print("noptimization finished, #iter = "+iter+"n");
  610. }
  612. // return 1 if already optimal, return 0 otherwise
  613. int select_working_set(int[] working_set)
  614. {
  615. // return i,j which maximize -grad(f)^T d , under constraint
  616. // if alpha_i == C, d != +1
  617. // if alpha_i == 0, d != -1
  619. double Gmax1 = -INF; // max { -grad(f)_i * d | y_i*d = +1 }
  620. int Gmax1_idx = -1;
  622. double Gmax2 = -INF; // max { -grad(f)_i * d | y_i*d = -1 }
  623. int Gmax2_idx = -1;
  625. for(int i=0;i<active_size;i++)
  626. {
  627. if(y[i]==+1) // y = +1
  628. {
  629. if(!is_upper_bound(i)) // d = +1
  630. {
  631. if(-G[i] > Gmax1)
  632. {
  633. Gmax1 = -G[i];
  634. Gmax1_idx = i;
  635. }
  636. }
  637. if(!is_lower_bound(i)) // d = -1
  638. {
  639. if(G[i] > Gmax2)
  640. {
  641. Gmax2 = G[i];
  642. Gmax2_idx = i;
  643. }
  644. }
  645. }
  646. else // y = -1
  647. {
  648. if(!is_upper_bound(i)) // d = +1
  649. {
  650. if(-G[i] > Gmax2)
  651. {
  652. Gmax2 = -G[i];
  653. Gmax2_idx = i;
  654. }
  655. }
  656. if(!is_lower_bound(i)) // d = -1
  657. {
  658. if(G[i] > Gmax1)
  659. {
  660. Gmax1 = G[i];
  661. Gmax1_idx = i;
  662. }
  663. }
  664. }
  665. }
  667. if(Gmax1+Gmax2 < eps)
  668.   return 1;
  670. working_set[0] = Gmax1_idx;
  671. working_set[1] = Gmax2_idx;
  672. return 0;
  673. }
  675. void do_shrinking()
  676. {
  677. int i,j,k;
  678. int[] working_set = new int[2];
  679. if(select_working_set(working_set)!=0) return;
  680. i = working_set[0];
  681. j = working_set[1];
  682. double Gm1 = -y[j]*G[j];
  683. double Gm2 = y[i]*G[i];
  685. // shrink
  687. for(k=0;k<active_size;k++)
  688. {
  689. if(is_lower_bound(k))
  690. {
  691. if(y[k]==+1)
  692. {
  693. if(-G[k] >= Gm1) continue;
  694. }
  695. else if(-G[k] >= Gm2) continue;
  696. }
  697. else if(is_upper_bound(k))
  698. {
  699. if(y[k]==+1)
  700. {
  701. if(G[k] >= Gm2) continue;
  702. }
  703. else if(G[k] >= Gm1) continue;
  704. }
  705. else continue;
  707. --active_size;
  708. swap_index(k,active_size);
  709. --k; // look at the newcomer
  710. }
  712. // unshrink, check all variables again before final iterations
  714. if(unshrinked || -(Gm1 + Gm2) > eps*10) return;
  716. unshrinked = true;
  717. reconstruct_gradient();
  719. for(k=l-1;k>=active_size;k--)
  720. {
  721. if(is_lower_bound(k))
  722. {
  723. if(y[k]==+1)
  724. {
  725. if(-G[k] < Gm1) continue;
  726. }
  727. else if(-G[k] < Gm2) continue;
  728. }
  729. else if(is_upper_bound(k))
  730. {
  731. if(y[k]==+1)
  732. {
  733. if(G[k] < Gm2) continue;
  734. }
  735. else if(G[k] < Gm1) continue;
  736. }
  737. else continue;
  739. swap_index(k,active_size);
  740. active_size++;
  741. ++k; // look at the newcomer
  742. }
  743. }
  745. double calculate_rho()
  746. {
  747. double r;
  748. int nr_free = 0;
  749. double ub = INF, lb = -INF, sum_free = 0;
  750. for(int i=0;i<active_size;i++)
  751. {
  752. double yG = y[i]*G[i];
  754. if(is_lower_bound(i))
  755. {
  756. if(y[i] > 0)
  757. ub = Math.min(ub,yG);
  758. else
  759. lb = Math.max(lb,yG);
  760. }
  761. else if(is_upper_bound(i))
  762. {
  763. if(y[i] < 0)
  764. ub = Math.min(ub,yG);
  765. else
  766. lb = Math.max(lb,yG);
  767. }
  768. else
  769. {
  770. ++nr_free;
  771. sum_free += yG;
  772. }
  773. }
  775. if(nr_free>0)
  776. r = sum_free/nr_free;
  777. else
  778. r = (ub+lb)/2;
  780. return r;
  781. }
  783. }
  785. //
  786. // Solver for nu-svm classification and regression
  787. //
  788. // additional constraint: e^T alpha = constant
  789. //
  790. final class Solver_NU extends Solver
  791. {
  792. private SolutionInfo si;
  794. void Solve(int l, Kernel Q, double[] b, byte[] y,
  795.    double[] alpha, double Cp, double Cn, double eps,
  796.    SolutionInfo si, int shrinking)
  797. {
  798. this.si = si;
  799. super.Solve(l,Q,b,y,alpha,Cp,Cn,eps,si,shrinking);
  800. }
  802. int select_working_set(int[] working_set)
  803. {
  804. // return i,j which maximize -grad(f)^T d , under constraint
  805. // if alpha_i == C, d != +1
  806. // if alpha_i == 0, d != -1
  808. double Gmax1 = -INF; // max { -grad(f)_i * d | y_i = +1, d = +1 }
  809. int Gmax1_idx = -1;
  811. double Gmax2 = -INF; // max { -grad(f)_i * d | y_i = +1, d = -1 }
  812. int Gmax2_idx = -1;
  814. double Gmax3 = -INF; // max { -grad(f)_i * d | y_i = -1, d = +1 }
  815. int Gmax3_idx = -1;
  817. double Gmax4 = -INF; // max { -grad(f)_i * d | y_i = -1, d = -1 }
  818. int Gmax4_idx = -1;
  820. for(int i=0;i<active_size;i++)
  821. {
  822. if(y[i]==+1) // y == +1
  823. {
  824. if(!is_upper_bound(i)) // d = +1
  825. {
  826. if(-G[i] > Gmax1)
  827. {
  828. Gmax1 = -G[i];
  829. Gmax1_idx = i;
  830. }
  831. }
  832. if(!is_lower_bound(i)) // d = -1
  833. {
  834. if(G[i] > Gmax2)
  835. {
  836. Gmax2 = G[i];
  837. Gmax2_idx = i;
  838. }
  839. }
  840. }
  841. else // y == -1
  842. {
  843. if(!is_upper_bound(i)) // d = +1
  844. {
  845. if(-G[i] > Gmax3)
  846. {
  847. Gmax3 = -G[i];
  848. Gmax3_idx = i;
  849. }
  850. }
  851. if(!is_lower_bound(i)) // d = -1
  852. {
  853. if(G[i] > Gmax4)
  854. {
  855. Gmax4 = G[i];
  856. Gmax4_idx = i;
  857. }
  858. }
  859. }
  860. }
  862. if(Math.max(Gmax1+Gmax2,Gmax3+Gmax4) < eps)
  863.   return 1;
  865. if(Gmax1+Gmax2 > Gmax3+Gmax4)
  866. {
  867. working_set[0] = Gmax1_idx;
  868. working_set[1] = Gmax2_idx;
  869. }
  870. else
  871. {
  872. working_set[0] = Gmax3_idx;
  873. working_set[1] = Gmax4_idx;
  874. }
  875. return 0;
  876. }
  878. void do_shrinking()
  879. {
  880. double Gmax1 = -INF; // max { -grad(f)_i * d | y_i = +1, d = +1 }
  881. double Gmax2 = -INF; // max { -grad(f)_i * d | y_i = +1, d = -1 }
  882. double Gmax3 = -INF; // max { -grad(f)_i * d | y_i = -1, d = +1 }
  883. double Gmax4 = -INF; // max { -grad(f)_i * d | y_i = -1, d = -1 }
  885. int k;
  886. for(k=0;k<active_size;k++)
  887. {
  888. if(!is_upper_bound(k))
  889. {
  890. if(y[k]==+1)
  891. {
  892. if(-G[k] > Gmax1) Gmax1 = -G[k];
  893. }
  894. else if(-G[k] > Gmax3) Gmax3 = -G[k];
  895. }
  896. if(!is_lower_bound(k))
  897. {
  898. if(y[k]==+1)
  899. {
  900. if(G[k] > Gmax2) Gmax2 = G[k];
  901. }
  902. else if(G[k] > Gmax4) Gmax4 = G[k];
  903. }
  904. }
  906. double Gm1 = -Gmax2;
  907. double Gm2 = -Gmax1;
  908. double Gm3 = -Gmax4;
  909. double Gm4 = -Gmax3;
  911. for(k=0;k<active_size;k++)
  912. {
  913. if(is_lower_bound(k))
  914. {
  915. if(y[k]==+1)
  916. {
  917. if(-G[k] >= Gm1) continue;
  918. }
  919. else if(-G[k] >= Gm3) continue;
  920. }
  921. else if(is_upper_bound(k))
  922. {
  923. if(y[k]==+1)
  924. {
  925. if(G[k] >= Gm2) continue;
  926. }
  927. else if(G[k] >= Gm4) continue;
  928. }
  929. else continue;
  931. --active_size;
  932. swap_index(k,active_size);
  933. --k; // look at the newcomer
  934. }
  936. // unshrink, check all variables again before final iterations
  938. if(unshrinked || Math.max(-(Gm1+Gm2),-(Gm3+Gm4)) > eps*10) return;
  940. unshrinked = true;
  941. reconstruct_gradient();
  943. for(k=l-1;k>=active_size;k--)
  944. {
  945. if(is_lower_bound(k))
  946. {
  947. if(y[k]==+1)
  948. {
  949. if(-G[k] < Gm1) continue;
  950. }
  951. else if(-G[k] < Gm3) continue;
  952. }
  953. else if(is_upper_bound(k))
  954. {
  955. if(y[k]==+1)
  956. {
  957. if(G[k] < Gm2) continue;
  958. }
  959. else if(G[k] < Gm4) continue;
  960. }
  961. else continue;
  963. swap_index(k,active_size);
  964. active_size++;
  965. ++k; // look at the newcomer
  966. }
  967. }
  969. double calculate_rho()
  970. {
  971. int nr_free1 = 0,nr_free2 = 0;
  972. double ub1 = INF, ub2 = INF;
  973. double lb1 = -INF, lb2 = -INF;
  974. double sum_free1 = 0, sum_free2 = 0;
  976. for(int i=0;i<active_size;i++)
  977. {
  978. if(y[i]==+1)
  979. {
  980. if(is_lower_bound(i))
  981. ub1 = Math.min(ub1,G[i]);
  982. else if(is_upper_bound(i))
  983. lb1 = Math.max(lb1,G[i]);
  984. else
  985. {
  986. ++nr_free1;
  987. sum_free1 += G[i];
  988. }
  989. }
  990. else
  991. {
  992. if(is_lower_bound(i))
  993. ub2 = Math.min(ub2,G[i]);
  994. else if(is_upper_bound(i))
  995. lb2 = Math.max(lb2,G[i]);
  996. else
  997. {
  998. ++nr_free2;
  999. sum_free2 += G[i];
  1000. }
  1001. }
  1002. }
  1004. double r1,r2;
  1005. if(nr_free1 > 0)
  1006. r1 = sum_free1/nr_free1;
  1007. else
  1008. r1 = (ub1+lb1)/2;
  1010. if(nr_free2 > 0)
  1011. r2 = sum_free2/nr_free2;
  1012. else
  1013. r2 = (ub2+lb2)/2;
  1015. si.r = (r1+r2)/2;
  1016. return (r1-r2)/2;
  1017. }
  1018. }
  1020. //
  1021. // Q matrices for various formulations
  1022. //
  1023. class SVC_Q extends Kernel
  1024. {
  1025. private final byte[] y;
  1026. private final Cache cache;
  1028. SVC_Q(svm_problem prob, svm_parameter param, byte[] y_)
  1029. {
  1030. super(prob.l, prob.x, param);
  1031. y = (byte[])y_.clone();
  1032. cache = new Cache(prob.l,(int)(param.cache_size*(1<<20)));
  1033. }
  1035. Qfloat[] get_Q(int i, int len)
  1036. {
  1037. Qfloat[][] data = new Qfloat[1][];
  1038. int start;
  1039. if((start = cache.get_data(i,data,len)) < len)
  1040. {
  1041. for(int j=start;j<len;j++)
  1042. data[0][j] = (Qfloat)(y[i]*y[j]*kernel_function(i,j));
  1043. }
  1044. return data[0];
  1045. }
  1047. void swap_index(int i, int j)
  1048. {
  1049. cache.swap_index(i,j);
  1050. super.swap_index(i,j);
  1051. swap(byte,y[i],y[j]);
  1052. }
  1053. }
  1055. class ONE_CLASS_Q extends Kernel
  1056. {
  1057. private final Cache cache;
  1059. ONE_CLASS_Q(svm_problem prob, svm_parameter param)
  1060. {
  1061. super(prob.l, prob.x, param);
  1062. cache = new Cache(prob.l,(int)(param.cache_size*(1<<20)));
  1063. }
  1065. Qfloat[] get_Q(int i, int len)
  1066. {
  1067. Qfloat[][] data = new Qfloat[1][];
  1068. int start;
  1069. if((start = cache.get_data(i,data,len)) < len)
  1070. {
  1071. for(int j=start;j<len;j++)
  1072. data[0][j] = (Qfloat)kernel_function(i,j);
  1073. }
  1074. return data[0];
  1075. }
  1077. void swap_index(int i, int j)
  1078. {
  1079. cache.swap_index(i,j);
  1080. super.swap_index(i,j);
  1081. }
  1082. }
  1084. class SVR_Q extends Kernel
  1085. {
  1086. private final int l;
  1087. private final Cache cache;
  1088. private final byte[] sign;
  1089. private final int[] index;
  1090. private int next_buffer;
  1091. private Qfloat[][] buffer;
  1093. SVR_Q(svm_problem prob, svm_parameter param)
  1094. {
  1095. super(prob.l, prob.x, param);
  1096. l = prob.l;
  1097. cache = new Cache(l,(int)(param.cache_size*(1<<20)));
  1098. sign = new byte[2*l];
  1099. index = new int[2*l];
  1100. for(int k=0;k<l;k++)
  1101. {
  1102. sign[k] = 1;
  1103. sign[k+l] = -1;
  1104. index[k] = k;
  1105. index[k+l] = k;
  1106. }
  1107. buffer = new Qfloat[2][2*l];
  1108. next_buffer = 0;
  1109. }
  1111. void swap_index(int i, int j)
  1112. {
  1113. swap(byte,sign[i],sign[j]);
  1114. swap(int,index[i],index[j]);
  1115. }
  1117. Qfloat[] get_Q(int i, int len)
  1118. {
  1119. Qfloat[][] data = new Qfloat[1][];
  1120. int real_i = index[i];
  1121. if(cache.get_data(real_i,data,l) < l)
  1122. {
  1123. for(int j=0;j<l;j++)
  1124. data[0][j] = (Qfloat)kernel_function(real_i,j);
  1125. }
  1127. // reorder and copy
  1128. Qfloat buf[] = buffer[next_buffer];
  1129. next_buffer = 1 - next_buffer;
  1130. byte si = sign[i];
  1131. for(int j=0;j<len;j++)
  1132. buf[j] = si * sign[j] * data[0][index[j]];
  1133. return buf;
  1134. }
  1135. }
  1137. public class svm {
  1138. //
  1139. // construct and solve various formulations
  1140. //
  1141. private static void solve_c_svc(svm_problem prob, svm_parameter param,
  1142. double[] alpha, Solver.SolutionInfo si,
  1143. double Cp, double Cn)
  1144. {
  1145. int l = prob.l;
  1146. double[] minus_ones = new double[l];
  1147. byte[] y = new byte[l];
  1149. int i;
  1151. for(i=0;i<l;i++)
  1152. {
  1153. alpha[i] = 0;
  1154. minus_ones[i] = -1;
  1155. if(prob.y[i] > 0) y[i] = +1; else y[i]=-1;
  1156. }
  1158. Solver s = new Solver();
  1159. s.Solve(l, new SVC_Q(prob,param,y), minus_ones, y,
  1160. alpha, Cp, Cn, param.eps, si, param.shrinking);
  1162. double sum_alpha=0;
  1163. for(i=0;i<l;i++)
  1164. sum_alpha += alpha[i];
  1166. System.out.print("nu = "+sum_alpha/(param.C*prob.l)+"n");
  1168. for(i=0;i<l;i++)
  1169. alpha[i] *= y[i];
  1170. }
  1172. private static void solve_nu_svc(svm_problem prob, svm_parameter param,
  1173.   double[] alpha, Solver.SolutionInfo si)
  1174. {
  1175. int i;
  1176. int l = prob.l;
  1177. double nu = param.nu;
  1179. int y_pos = 0;
  1180. int y_neg = 0;
  1181. byte[] y = new byte[l];
  1183. for(i=0;i<l;i++)
  1184. if(prob.y[i]>0)
  1185. {
  1186. y[i] = +1;
  1187. ++y_pos;
  1188. }
  1189. else
  1190. {
  1191. y[i] = -1;
  1192. ++y_neg;
  1193. }
  1195. if(nu < 0 || nu*l/2 > Math.min(y_pos,y_neg))
  1196. {
  1197. System.err.print("specified nu is infeasiblen");
  1198. System.exit(1);
  1199. }
  1201. double sum_pos = nu*l/2;
  1202. double sum_neg = nu*l/2;
  1204. for(i=0;i<l;i++)
  1205. if(y[i] == +1)
  1206. {
  1207. alpha[i] = Math.min(1.0,sum_pos);
  1208. sum_pos -= alpha[i];
  1209. }
  1210. else
  1211. {
  1212. alpha[i] = Math.min(1.0,sum_neg);
  1213. sum_neg -= alpha[i];
  1214. }
  1216. double[] zeros = new double[l];
  1218. for(i=0;i<l;i++)
  1219. zeros[i] = 0;
  1221. Solver_NU s = new Solver_NU();
  1222. s.Solve(l, new SVC_Q(prob,param,y), zeros, y,
  1223. alpha, 1.0, 1.0, param.eps, si, param.shrinking);
  1224. double r = si.r;
  1226. System.out.print("C = "+1/r+"n");
  1228. for(i=0;i<l;i++)
  1229. alpha[i] *= y[i]/r;
  1231. si.rho /= r;
  1232. si.obj /= (r*r);
  1233. si.upper_bound_p = 1/r;
  1234. si.upper_bound_n = 1/r;
  1235. }
  1237. private static void solve_one_class(svm_problem prob, svm_parameter param,
  1238.      double[] alpha, Solver.SolutionInfo si)
  1239. {
  1240. int l = prob.l;
  1241. double[] zeros = new double[l];
  1242. byte[] ones = new byte[l];
  1243. int i;
  1245. int n = (int)(param.nu*prob.l); // # of alpha's at upper bound
  1246. if(n>=prob.l)
  1247. {
  1248. System.err.print("nu must be in (0,1)n");
  1249. System.exit(1);
  1250. }
  1251. for(i=0;i<n;i++)
  1252. alpha[i] = 1;
  1253. alpha[n] = param.nu * prob.l - n;
  1254. for(i=n+1;i<l;i++)
  1255. alpha[i] = 0;
  1257. for(i=0;i<l;i++)
  1258. {
  1259. zeros[i] = 0;
  1260. ones[i] = 1;
  1261. }
  1263. Solver s = new Solver();
  1264. s.Solve(l, new ONE_CLASS_Q(prob,param), zeros, ones,
  1265. alpha, 1.0, 1.0, param.eps, si, param.shrinking);
  1266. }
  1268. private static void solve_epsilon_svr(svm_problem prob, svm_parameter param,
  1269. double[] alpha, Solver.SolutionInfo si)
  1270. {
  1271. int l = prob.l;
  1272. double[] alpha2 = new double[2*l];
  1273. double[] linear_term = new double[2*l];
  1274. byte[] y = new byte[2*l];
  1275. int i;
  1277. for(i=0;i<l;i++)
  1278. {
  1279. alpha2[i] = 0;
  1280. linear_term[i] = param.p - prob.y[i];
  1281. y[i] = 1;
  1283. alpha2[i+l] = 0;
  1284. linear_term[i+l] = param.p + prob.y[i];
  1285. y[i+l] = -1;
  1286. }
  1288. Solver s = new Solver();
  1289. s.Solve(2*l, new SVR_Q(prob,param), linear_term, y,
  1290. alpha2, param.C, param.C, param.eps, si, param.shrinking);
  1292. double sum_alpha = 0;
  1293. for(i=0;i<l;i++)
  1294. {
  1295. alpha[i] = alpha2[i] - alpha2[i+l];
  1296. sum_alpha += Math.abs(alpha[i]);
  1297. }
  1298. System.out.print("nu = "+sum_alpha/(param.C*l)+"n");
  1299. }
  1301. private static void solve_nu_svr(svm_problem prob, svm_parameter param,
  1302. double[] alpha, Solver.SolutionInfo si)
  1303. {
  1304. if(param.nu < 0 || param.nu > 1)
  1305. {
  1306. System.err.print("specified nu is out of rangen");
  1307. System.exit(1);
  1308. }
  1310. int l = prob.l;
  1311. double C = param.C;
  1312. double[] alpha2 = new double[2*l];
  1313. double[] linear_term = new double[2*l];
  1314. byte[] y = new byte[2*l];
  1315. int i;
  1317. double sum = C * param.nu * l / 2;
  1318. for(i=0;i<l;i++)
  1319. {
  1320. alpha2[i] = alpha2[i+l] = Math.min(sum,C);
  1321. sum -= alpha2[i];
  1323. linear_term[i] = - prob.y[i];
  1324. y[i] = 1;
  1326. linear_term[i+l] = prob.y[i];
  1327. y[i+l] = -1;
  1328. }
  1330. Solver_NU s = new Solver_NU();
  1331. s.Solve(2*l, new SVR_Q(prob,param), linear_term, y,
  1332. alpha2, C, C, param.eps, si, param.shrinking);
  1334. System.out.print("epsilon = "+(-si.r)+"n");
  1336. for(i=0;i<l;i++)
  1337. alpha[i] = alpha2[i] - alpha2[i+l];
  1338. }
  1340. //
  1341. // decision_function
  1342. //
  1343. static class decision_function
  1344. {
  1345. double[] alpha;
  1346. double rho;
  1347. };
  1349. static decision_function svm_train_one(
  1350. svm_problem prob, svm_parameter param,
  1351. double Cp, double Cn)
  1352. {
  1353. double[] alpha = new double[prob.l];
  1354. Solver.SolutionInfo si = new Solver.SolutionInfo();
  1355. switch(param.svm_type)
  1356. {
  1357. case svm_parameter.C_SVC:
  1358. solve_c_svc(prob,param,alpha,si,Cp,Cn);
  1359. break;
  1360. case svm_parameter.NU_SVC:
  1361. solve_nu_svc(prob,param,alpha,si);
  1362. break;
  1363. case svm_parameter.ONE_CLASS:
  1364. solve_one_class(prob,param,alpha,si);
  1365. break;
  1366. case svm_parameter.EPSILON_SVR:
  1367. solve_epsilon_svr(prob,param,alpha,si);
  1368. break;
  1369. case svm_parameter.NU_SVR:
  1370. solve_nu_svr(prob,param,alpha,si);
  1371. break;
  1372. }
  1374. System.out.print("obj = "+si.obj+", rho = "+si.rho+"n");
  1376. // output SVs
  1378. int nSV = 0;
  1379. int nBSV = 0;
  1380. for(int i=0;i<prob.l;i++)
  1381. {
  1382. if(Math.abs(alpha[i]) > 0)
  1383. {
  1384. ++nSV;
  1385. if(prob.y[i] > 0)
  1386. {
  1387. if(Math.abs(alpha[i]) >= si.upper_bound_p)
  1388. ++nBSV;
  1389. }
  1390. else
  1391. {
  1392. if(Math.abs(alpha[i]) >= si.upper_bound_n)
  1393. ++nBSV;
  1394. }
  1395. }
  1396. }
  1398. System.out.print("nSV = "+nSV+", nBSV = "+nBSV+"n");
  1400. decision_function f = new decision_function();
  1401. f.alpha = alpha;
  1402. f.rho = si.rho;
  1403. return f;
  1404. }
  1406. //
  1407. // Interface functions
  1408. //
  1409. public static svm_model svm_train(svm_problem prob, svm_parameter param)
  1410. {
  1411. svm_model model = new svm_model();
  1412. model.param = param;
  1414. if(param.svm_type == svm_parameter.ONE_CLASS ||
  1415.    param.svm_type == svm_parameter.EPSILON_SVR ||
  1416.    param.svm_type == svm_parameter.NU_SVR)
  1417. {
  1418. // regression or one-class-svm
  1419. model.nr_class = 2;
  1420. model.label = null;
  1421. model.nSV = null;
  1422. model.sv_coef = new double[1][];
  1423. decision_function f = svm_train_one(prob,param,0,0);
  1424. model.rho = new double[1];
  1425. model.rho[0] = f.rho;
  1427. int nSV = 0;
  1428. int i;
  1429. for(i=0;i<prob.l;i++)
  1430. if(Math.abs(f.alpha[i]) > 0) ++nSV;
  1431. model.l = nSV;
  1432. model.SV = new svm_node[nSV][];
  1433. model.sv_coef[0] = new double[nSV];
  1434. int j = 0;
  1435. for(i=0;i<prob.l;i++)
  1436. if(Math.abs(f.alpha[i]) > 0)
  1437. {
  1438. model.SV[j] = prob.x[i];
  1439. model.sv_coef[0][j] = f.alpha[i];
  1440. ++j;
  1441. }
  1442. }
  1443. else
  1444. {
  1445. // classification
  1446. // find out the number of classes
  1447. int l = prob.l;
  1448. int max_nr_class = 16;
  1449. int nr_class = 0;
  1450. int[] label = new int[max_nr_class];
  1451. int[] count = new int[max_nr_class];
  1452. int[] index = new int[l];
  1454. int i;
  1455. for(i=0;i<l;i++)
  1456. {
  1457. int this_label = (int)prob.y[i];
  1458. int j;
  1459. for(j=0;j<nr_class;j++)
  1460. if(this_label == label[j])
  1461. {
  1462. ++count[j];
  1463. break;
  1464. }
  1465. index[i] = j;
  1466. if(j == nr_class)
  1467. {
  1468. if(nr_class == max_nr_class)
  1469. {
  1470. max_nr_class *= 2;
  1471. int[] new_data = new int[max_nr_class];
  1472. System.arraycopy(label,0,new_data,0,label.length);
  1473. label = new_data;
  1475. new_data = new int[max_nr_class];
  1476. System.arraycopy(count,0,new_data,0,count.length);
  1477. count = new_data;
  1478. }
  1479. label[nr_class] = this_label;
  1480. count[nr_class] = 1;
  1481. ++nr_class;
  1482. }
  1483. }
  1485. // group training data of the same class
  1487. int[] start = new int[nr_class];
  1488. start[0] = 0;
  1489. for(i=1;i<nr_class;i++)
  1490. start[i] = start[i-1]+count[i-1];
  1492. svm_node[][] x = new svm_node[l][];
  1494. for(i=0;i<l;i++)
  1495. {
  1496. x[start[index[i]]] = prob.x[i];
  1497. ++start[index[i]];
  1498. }
  1500. start[0] = 0;
  1501. for(i=1;i<nr_class;i++)
  1502. start[i] = start[i-1]+count[i-1];
  1504. // calculate weighted C
  1506. double[] weighted_C = new double[nr_class];
  1507. for(i=0;i<nr_class;i++)
  1508. weighted_C[i] = param.C;
  1509. for(i=0;i<param.nr_weight;i++)
  1510. {
  1511. int j;
  1512. for(j=0;j<nr_class;j++)
  1513. if(param.weight_label[i] == label[j])
  1514. break;
  1515. if(j == nr_class)
  1516. System.err.print("warning: class label "+param.weight_label[i]+" specified in weight is not foundn");
  1517. else
  1518. weighted_C[j] *= param.weight[i];
  1519. }
  1521. // train n*(n-1)/2 models
  1523. boolean[] nonzero = new boolean[l];
  1524. for(i=0;i<l;i++)
  1525. nonzero[i] = false;
  1526. decision_function[] f = new decision_function[nr_class*(nr_class-1)/2];
  1528. int p = 0;
  1529. for(i=0;i<nr_class;i++)
  1530. for(int j=i+1;j<nr_class;j++)
  1531. {
  1532. svm_problem sub_prob = new svm_problem();
  1533. int si = start[i], sj = start[j];
  1534. int ci = count[i], cj = count[j];
  1535. sub_prob.l = ci+cj;
  1536. sub_prob.x = new svm_node[sub_prob.l][];
  1537. sub_prob.y = new double[sub_prob.l];
  1538. int k;
  1539. for(k=0;k<ci;k++)
  1540. {
  1541. sub_prob.x[k] = x[si+k];
  1542. sub_prob.y[k] = +1;
  1543. }
  1544. for(k=0;k<cj;k++)
  1545. {
  1546. sub_prob.x[ci+k] = x[sj+k];
  1547. sub_prob.y[ci+k] = -1;
  1548. }
  1550. f[p] = svm_train_one(sub_prob,param,weighted_C[i],weighted_C[j]);
  1551. for(k=0;k<ci;k++)
  1552. if(!nonzero[si+k] && Math.abs(f[p].alpha[k]) > 0)
  1553. nonzero[si+k] = true;
  1554. for(k=0;k<cj;k++)
  1555. if(!nonzero[sj+k] && Math.abs(f[p].alpha[ci+k]) > 0)
  1556. nonzero[sj+k] = true;
  1557. ++p;
  1558. }
  1560. // build output
  1562. model.nr_class = nr_class;
  1564. model.label = new int[nr_class];
  1565. for(i=0;i<nr_class;i++)
  1566. model.label[i] = label[i];
  1568. model.rho = new double[nr_class*(nr_class-1)/2];
  1569. for(i=0;i<nr_class*(nr_class-1)/2;i++)
  1570. model.rho[i] = f[i].rho;
  1572. int nnz = 0;
  1573. int[] nz_count = new int[nr_class];
  1574. model.nSV = new int[nr_class];
  1575. for(i=0;i<nr_class;i++)
  1576. {
  1577. int nSV = 0;
  1578. for(int j=0;j<count[i];j++)
  1579. if(nonzero[start[i]+j])
  1580. {
  1581. ++nSV;
  1582. ++nnz;
  1583. }
  1584. model.nSV[i] = nSV;
  1585. nz_count[i] = nSV;
  1586. }
  1588. System.out.print("Total nSV = "+nnz+"n");
  1590. model.l = nnz;
  1591. model.SV = new svm_node[nnz][];
  1592. p = 0;
  1593. for(i=0;i<l;i++)
  1594. if(nonzero[i]) model.SV[p++] = x[i];
  1596. int[] nz_start = new int[nr_class];
  1597. nz_start[0] = 0;
  1598. for(i=1;i<nr_class;i++)
  1599. nz_start[i] = nz_start[i-1]+nz_count[i-1];
  1601. model.sv_coef = new double[nr_class-1][];
  1602. for(i=0;i<nr_class-1;i++)
  1603. model.sv_coef[i] = new double[nnz];
  1605. p = 0;
  1606. for(i=0;i<nr_class;i++)
  1607. for(int j=i+1;j<nr_class;j++)
  1608. {
  1609. // classifier (i,j): coefficients with
  1610. // i are in sv_coef[j-1][nz_start[i]...],
  1611. // j are in sv_coef[i][nz_start[j]...]
  1613. int si = start[i];
  1614. int sj = start[j];
  1615. int ci = count[i];
  1616. int cj = count[j];
  1618. int q = nz_start[i];
  1619. int k;
  1620. for(k=0;k<ci;k++)
  1621. if(nonzero[si+k])
  1622. model.sv_coef[j-1][q++] = f[p].alpha[k];
  1623. q = nz_start[j];
  1624. for(k=0;k<cj;k++)
  1625. if(nonzero[sj+k])
  1626. model.sv_coef[i][q++] = f[p].alpha[ci+k];
  1627. ++p;
  1628. }
  1629. }
  1630. return model;
  1631. }
  1633. public static double svm_predict(svm_model model, svm_node[] x)
  1634. {
  1635. if(model.param.svm_type == svm_parameter.ONE_CLASS ||
  1636.    model.param.svm_type == svm_parameter.EPSILON_SVR ||
  1637.    model.param.svm_type == svm_parameter.NU_SVR)
  1638. {
  1639. double[] sv_coef = model.sv_coef[0];
  1640. double sum = 0;
  1641. for(int i=0;i<model.l;i++)
  1642. sum += sv_coef[i] * Kernel.k_function(x,model.SV[i],model.param);
  1643. sum -= model.rho[0];
  1644. if(model.param.svm_type == svm_parameter.ONE_CLASS)
  1645. return (sum>0)?1:-1;
  1646. else
  1647. return sum;
  1648. }
  1649. else
  1650. {
  1651. int i;
  1652. int nr_class = model.nr_class;
  1653. int l = model.l;
  1655. double[] kvalue = new double[l];
  1656. for(i=0;i<l;i++)
  1657. kvalue[i] = Kernel.k_function(x,model.SV[i],model.param);
  1659. int[] start = new int[nr_class];
  1660. start[0] = 0;
  1661. for(i=1;i<nr_class;i++)
  1662. start[i] = start[i-1]+model.nSV[i-1];
  1664. int[] vote = new int[nr_class];
  1665. for(i=0;i<nr_class;i++)
  1666. vote[i] = 0;
  1667. int p=0;
  1668. for(i=0;i<nr_class;i++)
  1669. for(int j=i+1;j<nr_class;j++)
  1670. {
  1671. double sum = 0;
  1672. int si = start[i];
  1673. int sj = start[j];
  1674. int ci = model.nSV[i];
  1675. int cj = model.nSV[j];
  1677. int k;
  1678. double[] coef1 = model.sv_coef[j-1];
  1679. double[] coef2 = model.sv_coef[i];
  1680. for(k=0;k<ci;k++)
  1681. sum += coef1[si+k] * kvalue[si+k];
  1682. for(k=0;k<cj;k++)
  1683. sum += coef2[sj+k] * kvalue[sj+k];
  1684. sum -= model.rho[p++];
  1685. if(sum > 0)
  1686. ++vote[i];
  1687. else
  1688. ++vote[j];
  1689. }
  1691. int vote_max_idx = 0;
  1692. for(i=1;i<nr_class;i++)
  1693. if(vote[i] > vote[vote_max_idx])
  1694. vote_max_idx = i;
  1695. return model.label[vote_max_idx];
  1696. }
  1697. }
  1698. static final String svm_type_table[] =
  1699. {
  1700. "c_svc","nu_svc","one_class","epsilon_svr","nu_svr",
  1701. };
  1703. static final String kernel_type_table[]=
  1704. {
  1705. "linear","polynomial","rbf","sigmoid",
  1706. };
  1708. public static void svm_save_model(String model_file_name, svm_model model) throws IOException
  1709. {
  1710. DataOutputStream fp = new DataOutputStream(new FileOutputStream(model_file_name));
  1712. svm_parameter param = model.param;
  1714. fp.writeBytes("svm_type "+svm_type_table[param.svm_type]+"n");
  1715. fp.writeBytes("kernel_type "+kernel_type_table[param.kernel_type]+"n");
  1717. if(param.kernel_type == svm_parameter.POLY)
  1718. fp.writeBytes("degree "+param.degree+"n");
  1720. if(param.kernel_type == svm_parameter.POLY ||
  1721.    param.kernel_type == svm_parameter.RBF ||
  1722.    param.kernel_type == svm_parameter.SIGMOID)
  1723. fp.writeBytes("gamma "+param.gamma+"n");
  1725. if(param.kernel_type == svm_parameter.POLY ||
  1726.    param.kernel_type == svm_parameter.SIGMOID)
  1727. fp.writeBytes("coef0 "+param.coef0+"n");
  1729. int nr_class = model.nr_class;
  1730. int l = model.l;
  1731. fp.writeBytes("nr_class "+nr_class+"n");
  1732. fp.writeBytes("total_sv "+l+"n");
  1734. {
  1735. fp.writeBytes("rho");
  1736. for(int i=0;i<nr_class*(nr_class-1)/2;i++)
  1737. fp.writeBytes(" "+model.rho[i]);
  1738. fp.writeBytes("n");
  1739. }
  1741. if(model.label != null)
  1742. {
  1743. fp.writeBytes("label");
  1744. for(int i=0;i<nr_class;i++)
  1745. fp.writeBytes(" "+model.label[i]);
  1746. fp.writeBytes("n");
  1747. }
  1749. if(model.nSV != null)
  1750. {
  1751. fp.writeBytes("nr_sv");
  1752. for(int i=0;i<nr_class;i++)
  1753. fp.writeBytes(" "+model.nSV[i]);
  1754. fp.writeBytes("n");
  1755. }
  1757. fp.writeBytes("SVn");
  1758. double[][] sv_coef = model.sv_coef;
  1759. svm_node[][] SV = model.SV;
  1761. for(int i=0;i<l;i++)
  1762. {
  1763. for(int j=0;j<nr_class-1;j++)
  1764. fp.writeBytes(sv_coef[j][i]+" ");
  1766. svm_node[] p = SV[i];
  1767. for(int j=0;j<p.length;j++)
  1768. fp.writeBytes(p[j].index+":"+p[j].value+" ");
  1769. fp.writeBytes("n");
  1770. }
  1772. fp.close();
  1773. }
  1775. private static double atof(String s)
  1776. {
  1777. return Double.valueOf(s).doubleValue();
  1778. }
  1780. private static int atoi(String s)
  1781. {
  1782. return Integer.parseInt(s);
  1783. }
  1785. public static svm_model svm_load_model(String model_file_name) throws IOException
  1786. {
  1787. BufferedReader fp = new BufferedReader(new FileReader(model_file_name));
  1789. // read parameters
  1791. svm_model model = new svm_model();
  1792. svm_parameter param = new svm_parameter();
  1793. model.param = param;
  1794. model.label = null;
  1795. model.nSV = null;
  1797. while(true)
  1798. {
  1799. String cmd = fp.readLine();
  1800. String arg = cmd.substring(cmd.indexOf(' ')+1);
  1802. if(cmd.startsWith("svm_type"))
  1803. {
  1804. int i;
  1805. for(i=0;i<svm_type_table.length;i++)
  1806. {
  1807. if(arg.indexOf(svm_type_table[i])!=-1)
  1808. {
  1809. param.svm_type=i;
  1810. break;
  1811. }
  1812. }
  1813. if(i == svm_type_table.length)
  1814. {
  1815. System.err.print("unknown svm type.n");
  1816. System.exit(1);
  1817. }
  1818. }
  1819. else if(cmd.startsWith("kernel_type"))
  1820. {
  1821. int i;
  1822. for(i=0;i<kernel_type_table.length;i++)
  1823. {
  1824. if(arg.indexOf(kernel_type_table[i])!=-1)
  1825. {
  1826. param.kernel_type=i;
  1827. break;
  1828. }
  1829. }
  1830. if(i == kernel_type_table.length)
  1831. {
  1832. System.err.print("unknown kernel function.n");
  1833. System.exit(1);
  1834. }
  1835. }
  1836. else if(cmd.startsWith("degree"))
  1837. param.degree = atof(arg);
  1838. else if(cmd.startsWith("gamma"))
  1839. param.gamma = atof(arg);
  1840. else if(cmd.startsWith("coef0"))
  1841. param.coef0 = atof(arg);
  1842. else if(cmd.startsWith("nr_class"))
  1843. model.nr_class = atoi(arg);
  1844. else if(cmd.startsWith("total_sv"))
  1845. model.l = atoi(arg);
  1846. else if(cmd.startsWith("rho"))
  1847. {
  1848. int n = model.nr_class * (model.nr_class-1)/2;
  1849. model.rho = new double[n];
  1850. StringTokenizer st = new StringTokenizer(arg);
  1851. for(int i=0;i<n;i++)
  1852. model.rho[i] = atof(st.nextToken());
  1853. }
  1854. else if(cmd.startsWith("label"))
  1855. {
  1856. int n = model.nr_class;
  1857. model.label = new int[n];
  1858. StringTokenizer st = new StringTokenizer(arg);
  1859. for(int i=0;i<n;i++)
  1860. model.label[i] = atoi(st.nextToken());
  1861. }
  1862. else if(cmd.startsWith("nr_sv"))
  1863. {
  1864. int n = model.nr_class;
  1865. model.nSV = new int[n];
  1866. StringTokenizer st = new StringTokenizer(arg);
  1867. for(int i=0;i<n;i++)
  1868. model.nSV[i] = atoi(st.nextToken());
  1869. }
  1870. else if(cmd.startsWith("SV"))
  1871. {
  1872. break;
  1873. }
  1874. else
  1875. {
  1876. System.err.print("unknown text in model filen");
  1877. System.exit(1);
  1878. }
  1879. }
  1881. // read sv_coef and SV
  1883. int m = model.nr_class - 1;
  1884. int l = model.l;
  1885. model.sv_coef = new double[m][l];
  1886. model.SV = new svm_node[l][];
  1888. for(int i=0;i<l;i++)
  1889. {
  1890. String line = fp.readLine();
  1891. StringTokenizer st = new StringTokenizer(line," tnrf:");
  1893. for(int k=0;k<m;k++)
  1894. model.sv_coef[k][i] = atof(st.nextToken());
  1895. int n = st.countTokens()/2;
  1896. model.SV[i] = new svm_node[n];
  1897. for(int j=0;j<n;j++)
  1898. {
  1899. model.SV[i][j] = new svm_node();
  1900. model.SV[i][j].index = atoi(st.nextToken());
  1901. model.SV[i][j].value = atof(st.nextToken());
  1902. }
  1903. }
  1905. fp.close();
  1906. return model;
  1907. }
  1908. }