Changeset 2415 for trunk/sources/LibSVM/Solver.cs
- Timestamp:
- 10/07/09 11:58:21 (15 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/sources/LibSVM/Solver.cs
r1819 r2415 19 19 20 20 using System; 21 using System.Linq; 21 22 using System.Collections.Generic; 22 23 using System.Diagnostics; 24 using System.IO; 23 25 24 26 namespace SVM 25 27 { 26 //27 // Kernel evaluation28 //29 // the static method k_function is for doing single kernel evaluation30 // the constructor of Kernel prepares to calculate the l*l kernel matrix31 // the member function get_Q is for getting one column from the Q Matrix32 //33 internal abstract class QMatrix34 {35 public abstract float[] get_Q(int column, int len);36 public abstract float[] get_QD();37 public abstract void swap_index(int i, int j);38 }39 40 internal abstract class Kernel : QMatrix41 {42 private Node[][] _x;43 private double[] _x_square;44 45 // Parameter46 private KernelType kernel_type;47 private int degree;48 private double gamma;49 private double coef0;50 51 public override void swap_index(int i, int j)52 {53 do { Node[] _ = _x[i]; _x[i] = _x[j]; _x[j] = _; } while (false);54 if (_x_square != null) do { double _ = _x_square[i]; _x_square[i] = _x_square[j]; _x_square[j] = _; } while (false);55 }56 57 private static double powi(double baseValue, int times)58 {59 double tmp = baseValue, ret = 1.0;60 61 for (int t = times; t > 0; t /= 2)62 {63 if (t % 2 == 1) ret *= tmp;64 tmp = tmp * tmp;65 }66 return ret;67 }68 69 private static double tanh(double x)70 {71 double e = Math.Exp(x);72 return 1.0 - 2.0 / (e * e + 1);73 }74 75 public double kernel_function(int i, int j)76 {77 switch (kernel_type)78 {79 case KernelType.LINEAR:80 return dot(_x[i], _x[j]);81 case KernelType.POLY:82 return powi(gamma * dot(_x[i], _x[j]) + coef0, degree);83 case KernelType.RBF:84 return Math.Exp(-gamma * (_x_square[i] + _x_square[j] - 2 * dot(_x[i], _x[j])));85 case KernelType.SIGMOID:86 return tanh(gamma * dot(_x[i], _x[j]) + coef0);87 case KernelType.PRECOMPUTED:88 return _x[i][(int)(_x[j][0].Value)].Value;89 default:90 return 0;91 }92 }93 94 public Kernel(int l, Node[][] x_, Parameter param)95 {96 this.kernel_type = param.KernelType;97 this.degree = param.Degree;98 this.gamma = param.Gamma;99 this.coef0 = param.Coefficient0;100 101 _x = (Node[][])x_.Clone();102 103 if (kernel_type == KernelType.RBF)104 {105 _x_square = new double[l];106 for (int i = 0; i < l; i++)107 _x_square[i] = dot(_x[i], _x[i]);108 }109 else _x_square = null;110 }111 112 public static double dot(Node[] x, Node[] y)113 {114 double sum = 0;115 int xlen = x.Length;116 int ylen = y.Length;117 int i = 0;118 int j = 0;119 while (i < xlen && j < ylen)120 {121 if (x[i].Index == y[j].Index)122 sum += x[i++].Value * y[j++].Value;123 else124 {125 if (x[i].Index > y[j].Index)126 ++j;127 else128 ++i;129 }130 }131 return sum;132 }133 134 public static double k_function(Node[] x, Node[] y, Parameter param)135 {136 switch (param.KernelType)137 {138 case KernelType.LINEAR:139 return dot(x, y);140 case KernelType.POLY:141 return powi(param.Gamma * dot(x, y) + param.Coefficient0, param.Degree);142 case KernelType.RBF:143 {144 double sum = 0;145 int xlen = x.Length;146 int ylen = y.Length;147 int i = 0;148 int j = 0;149 while (i < xlen && j < ylen)150 {151 if (x[i].Index == y[j].Index)152 {153 double d = x[i++].Value - y[j++].Value;154 sum += d * d;155 }156 else if (x[i].Index > y[j].Index)157 {158 sum += y[j].Value * y[j].Value;159 ++j;160 }161 else162 {163 sum += x[i].Value * x[i].Value;164 ++i;165 }166 }167 168 while (i < xlen)169 {170 sum += x[i].Value * x[i].Value;171 ++i;172 }173 174 while (j < ylen)175 {176 sum += y[j].Value * y[j].Value;177 ++j;178 }179 180 return Math.Exp(-param.Gamma * sum);181 }182 case KernelType.SIGMOID:183 return tanh(param.Gamma * dot(x, y) + param.Coefficient0);184 case KernelType.PRECOMPUTED:185 return x[(int)(y[0].Value)].Value;186 default:187 return 0;188 }189 }190 }191 28 192 29 // An SMO algorithm in Fan et al., JMLR 6(2005), p. 1889--1918 193 30 // Solves: 194 31 // 195 // min 0.5(\alpha^T Q \alpha) + p^T \alpha32 // Min 0.5(\alpha^T Q \alpha) + p^T \alpha 196 33 // 197 34 // y^T \alpha = \delta … … 211 48 { 212 49 protected int active_size; 213 protected s hort[] y;50 protected sbyte[] y; 214 51 protected double[] G; // gradient of objective function 215 pr otectedconst byte LOWER_BOUND = 0;216 pr otectedconst byte UPPER_BOUND = 1;217 pr otectedconst byte FREE = 2;218 pr otectedbyte[] alpha_status; // LOWER_BOUND, UPPER_BOUND, FREE219 pr otecteddouble[] alpha;220 protected QMatrix Q;52 private const byte LOWER_BOUND = 0; 53 private const byte UPPER_BOUND = 1; 54 private const byte FREE = 2; 55 private byte[] alpha_status; // LOWER_BOUND, UPPER_BOUND, FREE 56 private double[] alpha; 57 protected IQMatrix Q; 221 58 protected float[] QD; 222 protected double eps;223 pr otecteddouble Cp, Cn;224 pr otecteddouble[] p;225 pr otectedint[] active_set;226 pr otecteddouble[] G_bar; // gradient, if we treat free variables as 059 protected double EPS; 60 private double Cp, Cn; 61 private double[] p; 62 private int[] active_set; 63 private double[] G_bar; // gradient, if we treat free variables as 0 227 64 protected int l; 228 protected bool unshrink ed; // XXX65 protected bool unshrink; // XXX 229 66 230 67 protected const double INF = double.PositiveInfinity; 231 68 232 pr otecteddouble get_C(int i)69 private double get_C(int i) 233 70 { 234 71 return (y[i] > 0) ? Cp : Cn; 235 72 } 236 protected void update_alpha_status(int i) 73 74 private void update_alpha_status(int i) 237 75 { 238 76 if (alpha[i] >= get_C(i)) … … 242 80 else alpha_status[i] = FREE; 243 81 } 82 244 83 protected bool is_upper_bound(int i) { return alpha_status[i] == UPPER_BOUND; } 245 84 protected bool is_lower_bound(int i) { return alpha_status[i] == LOWER_BOUND; } 246 protected bool is_free(int i) { return alpha_status[i] == FREE; } 247 248 // java: information about solution except alpha, 249 // because we cannot return multiple values otherwise... 250 internal class SolutionInfo 85 86 private bool is_free(int i) { return alpha_status[i] == FREE; } 87 88 public class SolutionInfo 251 89 { 252 90 public double obj; … … 259 97 protected void swap_index(int i, int j) 260 98 { 261 Q. swap_index(i, j);262 do { short _ = y[i]; y[i] = y[j]; y[j] = _; } while (false);263 do { double _ = G[i]; G[i] = G[j]; G[j] = _; } while (false);264 do { byte _ = alpha_status[i]; alpha_status[i] = alpha_status[j]; alpha_status[j] = _; } while (false);265 do { double _ = alpha[i]; alpha[i] = alpha[j]; alpha[j] = _; } while (false);266 do { double _ = p[i]; p[i] = p[j]; p[j] = _; } while (false);267 do { int _ = active_set[i]; active_set[i] = active_set[j]; active_set[j] = _; } while (false);268 do { double _ = G_bar[i]; G_bar[i] = G_bar[j]; G_bar[j] = _; } while (false);99 Q.SwapIndex(i, j); 100 y.SwapIndex(i, j); 101 G.SwapIndex(i, j); 102 alpha_status.SwapIndex(i, j); 103 alpha.SwapIndex(i, j); 104 p.SwapIndex(i, j); 105 active_set.SwapIndex(i, j); 106 G_bar.SwapIndex(i, j); 269 107 } 270 108 … … 275 113 if (active_size == l) return; 276 114 277 int i; 278 for (i = active_size; i < l; i++) 279 G[i] = G_bar[i] + p[i]; 280 281 for (i = 0; i < active_size; i++) 282 if (is_free(i)) 283 { 284 float[] Q_i = Q.get_Q(i, l); 285 double alpha_i = alpha[i]; 286 for (int j = active_size; j < l; j++) 287 G[j] += alpha_i * Q_i[j]; 288 } 289 } 290 291 public virtual void Solve(int l, QMatrix Q, double[] p_, short[] y_, 292 double[] alpha_, double Cp, double Cn, double eps, SolutionInfo si, bool shrinking) 115 int i, j; 116 int nr_free = 0; 117 118 for (j = active_size; j < l; j++) 119 G[j] = G_bar[j] + p[j]; 120 121 for (j = 0; j < active_size; j++) 122 if (is_free(j)) 123 nr_free++; 124 125 if (2 * nr_free < active_size) 126 Procedures.info("\nWarning: using -h 0 may be faster\n"); 127 128 if (nr_free * l > 2 * active_size * (l - active_size)) 129 { 130 for (i = active_size; i < l; i++) 131 { 132 float[] Q_i = Q.GetQ(i, active_size); 133 for (j = 0; j < active_size; j++) 134 if (is_free(j)) 135 G[i] += alpha[j] * Q_i[j]; 136 } 137 } 138 else 139 { 140 for (i = 0; i < active_size; i++) 141 if (is_free(i)) 142 { 143 float[] Q_i = Q.GetQ(i, l); 144 double alpha_i = alpha[i]; 145 for (j = active_size; j < l; j++) 146 G[j] += alpha_i * Q_i[j]; 147 } 148 } 149 } 150 151 public virtual void Solve(int l, IQMatrix Q, double[] p_, sbyte[] y_, double[] alpha_, double Cp, double Cn, double eps, SolutionInfo si, bool shrinking) 293 152 { 294 153 this.l = l; 295 154 this.Q = Q; 296 QD = Q. get_QD();155 QD = Q.GetQD(); 297 156 p = (double[])p_.Clone(); 298 y = (s hort[])y_.Clone();157 y = (sbyte[])y_.Clone(); 299 158 alpha = (double[])alpha_.Clone(); 300 159 this.Cp = Cp; 301 160 this.Cn = Cn; 302 this. eps= eps;303 this.unshrink ed= false;161 this.EPS = eps; 162 this.unshrink = false; 304 163 305 164 // initialize alpha_status … … 331 190 if (!is_lower_bound(i)) 332 191 { 333 float[] Q_i = Q. get_Q(i, l);192 float[] Q_i = Q.GetQ(i, l); 334 193 double alpha_i = alpha[i]; 335 194 int j; … … 356 215 counter = Math.Min(l, 1000); 357 216 if (shrinking) do_shrinking(); 358 Debug.Write(".");217 Procedures.info("."); 359 218 } 360 219 … … 365 224 // reset active set size and check 366 225 active_size = l; 367 Debug.Write("*");226 Procedures.info("*"); 368 227 if (select_working_set(working_set) != 0) 369 228 break; … … 379 238 // update alpha[i] and alpha[j], handle bounds carefully 380 239 381 float[] Q_i = Q. get_Q(i, active_size);382 float[] Q_j = Q. get_Q(j, active_size);240 float[] Q_i = Q.GetQ(i, active_size); 241 float[] Q_j = Q.GetQ(j, active_size); 383 242 384 243 double C_i = get_C(i); … … 495 354 if (ui != is_upper_bound(i)) 496 355 { 497 Q_i = Q. get_Q(i, l);356 Q_i = Q.GetQ(i, l); 498 357 if (ui) 499 358 for (k = 0; k < l; k++) … … 506 365 if (uj != is_upper_bound(j)) 507 366 { 508 Q_j = Q. get_Q(j, l);367 Q_j = Q.GetQ(j, l); 509 368 if (uj) 510 369 for (k = 0; k < l; k++) … … 541 400 si.upper_bound_n = Cn; 542 401 543 Debug.Write("\noptimization finished, #iter = " + iter + "\n");402 Procedures.info("\noptimization finished, #iter = " + iter + "\n"); 544 403 } 545 404 546 405 // return 1 if already optimal, return 0 otherwise 547 protected virtualint select_working_set(int[] working_set)406 int select_working_set(int[] working_set) 548 407 { 549 408 // return i,j such that 550 // i: maximizes -y_i * grad(f)_i, i in I_up(\alpha)409 // i: Maximizes -y_i * grad(f)_i, i in I_up(\alpha) 551 410 // j: mimimizes the decrease of obj value 552 411 // (if quadratic coefficeint <= 0, replace it with tau) 553 412 // -y_j*grad(f)_j < -y_i*grad(f)_i, j in I_low(\alpha) 554 413 555 double G max = -INF;556 double G max2 = -INF;557 int G max_idx = -1;558 int G min_idx = -1;559 double obj_diff_ min = INF;414 double GMax = -INF; 415 double GMax2 = -INF; 416 int GMax_idx = -1; 417 int GMin_idx = -1; 418 double obj_diff_Min = INF; 560 419 561 420 for (int t = 0; t < active_size; t++) … … 563 422 { 564 423 if (!is_upper_bound(t)) 565 if (-G[t] >= G max)566 { 567 G max = -G[t];568 G max_idx = t;424 if (-G[t] >= GMax) 425 { 426 GMax = -G[t]; 427 GMax_idx = t; 569 428 } 570 429 } … … 572 431 { 573 432 if (!is_lower_bound(t)) 574 if (G[t] >= G max)575 { 576 G max = G[t];577 G max_idx = t;578 } 579 } 580 581 int i = G max_idx;433 if (G[t] >= GMax) 434 { 435 GMax = G[t]; 436 GMax_idx = t; 437 } 438 } 439 440 int i = GMax_idx; 582 441 float[] Q_i = null; 583 if (i != -1) // null Q_i not accessed: G max=-INF if i=-1584 Q_i = Q. get_Q(i, active_size);442 if (i != -1) // null Q_i not accessed: GMax=-INF if i=-1 443 Q_i = Q.GetQ(i, active_size); 585 444 586 445 for (int j = 0; j < active_size; j++) … … 590 449 if (!is_lower_bound(j)) 591 450 { 592 double grad_diff = G max + G[j];593 if (G[j] >= G max2)594 G max2 = G[j];451 double grad_diff = GMax + G[j]; 452 if (G[j] >= GMax2) 453 GMax2 = G[j]; 595 454 if (grad_diff > 0) 596 455 { 597 456 double obj_diff; 598 double quad_coef = Q_i[i] + QD[j] - 2 * y[i] * Q_i[j];457 double quad_coef = Q_i[i] + QD[j] - 2.0 * y[i] * Q_i[j]; 599 458 if (quad_coef > 0) 600 459 obj_diff = -(grad_diff * grad_diff) / quad_coef; … … 602 461 obj_diff = -(grad_diff * grad_diff) / 1e-12; 603 462 604 if (obj_diff <= obj_diff_ min)463 if (obj_diff <= obj_diff_Min) 605 464 { 606 G min_idx = j;607 obj_diff_ min = obj_diff;465 GMin_idx = j; 466 obj_diff_Min = obj_diff; 608 467 } 609 468 } … … 614 473 if (!is_upper_bound(j)) 615 474 { 616 double grad_diff = G max - G[j];617 if (-G[j] >= G max2)618 G max2 = -G[j];475 double grad_diff = GMax - G[j]; 476 if (-G[j] >= GMax2) 477 GMax2 = -G[j]; 619 478 if (grad_diff > 0) 620 479 { 621 480 double obj_diff; 622 double quad_coef = Q_i[i] + QD[j] + 2 * y[i] * Q_i[j];481 double quad_coef = Q_i[i] + QD[j] + 2.0 * y[i] * Q_i[j]; 623 482 if (quad_coef > 0) 624 483 obj_diff = -(grad_diff * grad_diff) / quad_coef; … … 626 485 obj_diff = -(grad_diff * grad_diff) / 1e-12; 627 486 628 if (obj_diff <= obj_diff_ min)487 if (obj_diff <= obj_diff_Min) 629 488 { 630 G min_idx = j;631 obj_diff_ min = obj_diff;489 GMin_idx = j; 490 obj_diff_Min = obj_diff; 632 491 } 633 492 } … … 636 495 } 637 496 638 if (G max + Gmax2 < eps)497 if (GMax + GMax2 < EPS) 639 498 return 1; 640 499 641 working_set[0] = G max_idx;642 working_set[1] = G min_idx;500 working_set[0] = GMax_idx; 501 working_set[1] = GMin_idx; 643 502 return 0; 644 503 } 645 504 646 private bool be_shrunk en(int i, double Gmax1, double Gmax2)505 private bool be_shrunk(int i, double GMax1, double GMax2) 647 506 { 648 507 if (is_upper_bound(i)) 649 508 { 650 509 if (y[i] == +1) 651 return (-G[i] > G max1);510 return (-G[i] > GMax1); 652 511 else 653 return (-G[i] > G max2);512 return (-G[i] > GMax2); 654 513 } 655 514 else if (is_lower_bound(i)) 656 515 { 657 516 if (y[i] == +1) 658 return (G[i] > G max2);517 return (G[i] > GMax2); 659 518 else 660 return (G[i] > G max1);519 return (G[i] > GMax1); 661 520 } 662 521 else … … 664 523 } 665 524 666 protected virtualvoid do_shrinking()525 void do_shrinking() 667 526 { 668 527 int i; 669 double G max1 = -INF; // max { -y_i * grad(f)_i | i in I_up(\alpha) }670 double G max2 = -INF; // max { y_i * grad(f)_i | i in I_low(\alpha) }671 672 // find maximal violating pair first528 double GMax1 = -INF; // Max { -y_i * grad(f)_i | i in I_up(\alpha) } 529 double GMax2 = -INF; // Max { y_i * grad(f)_i | i in I_low(\alpha) } 530 531 // find Maximal violating pair first 673 532 for (i = 0; i < active_size; i++) 674 533 { … … 677 536 if (!is_upper_bound(i)) 678 537 { 679 if (-G[i] >= G max1)680 G max1 = -G[i];538 if (-G[i] >= GMax1) 539 GMax1 = -G[i]; 681 540 } 682 541 if (!is_lower_bound(i)) 683 542 { 684 if (G[i] >= G max2)685 G max2 = G[i];543 if (G[i] >= GMax2) 544 GMax2 = G[i]; 686 545 } 687 546 } … … 690 549 if (!is_upper_bound(i)) 691 550 { 692 if (-G[i] >= G max2)693 G max2 = -G[i];551 if (-G[i] >= GMax2) 552 GMax2 = -G[i]; 694 553 } 695 554 if (!is_lower_bound(i)) 696 555 { 697 if (G[i] >= Gmax1) 698 Gmax1 = G[i]; 699 } 700 } 701 } 702 703 // shrink 556 if (G[i] >= GMax1) 557 GMax1 = G[i]; 558 } 559 } 560 } 561 562 if (unshrink == false && GMax1 + GMax2 <= EPS * 10) 563 { 564 unshrink = true; 565 reconstruct_gradient(); 566 active_size = l; 567 } 704 568 705 569 for (i = 0; i < active_size; i++) 706 if (be_shrunk en(i, Gmax1, Gmax2))570 if (be_shrunk(i, GMax1, GMax2)) 707 571 { 708 572 active_size--; 709 573 while (active_size > i) 710 574 { 711 if (!be_shrunk en(active_size, Gmax1, Gmax2))575 if (!be_shrunk(active_size, GMax1, GMax2)) 712 576 { 713 577 swap_index(i, active_size); … … 717 581 } 718 582 } 719 720 // unshrink, check all variables again before sealed iterations 721 722 if (unshrinked || Gmax1 + Gmax2 > eps * 10) return; 723 724 unshrinked = true; 725 reconstruct_gradient(); 726 727 for (i = l - 1; i >= active_size; i--) 728 if (!be_shrunken(i, Gmax1, Gmax2)) 729 { 730 while (active_size < i) 731 { 732 if (be_shrunken(active_size, Gmax1, Gmax2)) 733 { 734 swap_index(i, active_size); 735 break; 736 } 737 active_size++; 738 } 739 active_size++; 740 } 741 } 742 743 protected virtual double calculate_rho() 583 } 584 585 double calculate_rho() 744 586 { 745 587 double r; … … 786 628 // additional constraint: e^T \alpha = constant 787 629 // 788 sealedclass Solver_NU : Solver630 class Solver_NU : Solver 789 631 { 790 632 private SolutionInfo si; 791 633 792 public override void Solve(int l, QMatrix Q, double[] p, short[] y,634 public sealed override void Solve(int l, IQMatrix Q, double[] p, sbyte[] y, 793 635 double[] alpha, double Cp, double Cn, double eps, 794 636 SolutionInfo si, bool shrinking) … … 799 641 800 642 // return 1 if already optimal, return 0 otherwise 801 pr otected override int select_working_set(int[] working_set)643 private int select_working_set(int[] working_set) 802 644 { 803 645 // return i,j such that y_i = y_j and 804 // i: maximizes -y_i * grad(f)_i, i in I_up(\alpha)805 // j: minimizes the decrease of obj value646 // i: Maximizes -y_i * grad(f)_i, i in I_up(\alpha) 647 // j: Minimizes the decrease of obj value 806 648 // (if quadratic coefficeint <= 0, replace it with tau) 807 649 // -y_j*grad(f)_j < -y_i*grad(f)_i, j in I_low(\alpha) 808 650 809 double G maxp = -INF;810 double G maxp2 = -INF;811 int G maxp_idx = -1;812 813 double G maxn = -INF;814 double G maxn2 = -INF;815 int G maxn_idx = -1;816 817 int G min_idx = -1;818 double obj_diff_ min = INF;651 double GMaxp = -INF; 652 double GMaxp2 = -INF; 653 int GMaxp_idx = -1; 654 655 double GMaxn = -INF; 656 double GMaxn2 = -INF; 657 int GMaxn_idx = -1; 658 659 int GMin_idx = -1; 660 double obj_diff_Min = INF; 819 661 820 662 for (int t = 0; t < active_size; t++) … … 822 664 { 823 665 if (!is_upper_bound(t)) 824 if (-G[t] >= G maxp)825 { 826 G maxp = -G[t];827 G maxp_idx = t;666 if (-G[t] >= GMaxp) 667 { 668 GMaxp = -G[t]; 669 GMaxp_idx = t; 828 670 } 829 671 } … … 831 673 { 832 674 if (!is_lower_bound(t)) 833 if (G[t] >= G maxn)834 { 835 G maxn = G[t];836 G maxn_idx = t;837 } 838 } 839 840 int ip = G maxp_idx;841 int iN = G maxn_idx;675 if (G[t] >= GMaxn) 676 { 677 GMaxn = G[t]; 678 GMaxn_idx = t; 679 } 680 } 681 682 int ip = GMaxp_idx; 683 int iN = GMaxn_idx; 842 684 float[] Q_ip = null; 843 685 float[] Q_in = null; 844 if (ip != -1) // null Q_ip not accessed: G maxp=-INF if ip=-1845 Q_ip = Q. get_Q(ip, active_size);686 if (ip != -1) // null Q_ip not accessed: GMaxp=-INF if ip=-1 687 Q_ip = Q.GetQ(ip, active_size); 846 688 if (iN != -1) 847 Q_in = Q. get_Q(iN, active_size);689 Q_in = Q.GetQ(iN, active_size); 848 690 849 691 for (int j = 0; j < active_size; j++) … … 853 695 if (!is_lower_bound(j)) 854 696 { 855 double grad_diff = G maxp + G[j];856 if (G[j] >= G maxp2)857 G maxp2 = G[j];697 double grad_diff = GMaxp + G[j]; 698 if (G[j] >= GMaxp2) 699 GMaxp2 = G[j]; 858 700 if (grad_diff > 0) 859 701 { … … 865 707 obj_diff = -(grad_diff * grad_diff) / 1e-12; 866 708 867 if (obj_diff <= obj_diff_ min)709 if (obj_diff <= obj_diff_Min) 868 710 { 869 G min_idx = j;870 obj_diff_ min = obj_diff;711 GMin_idx = j; 712 obj_diff_Min = obj_diff; 871 713 } 872 714 } … … 877 719 if (!is_upper_bound(j)) 878 720 { 879 double grad_diff = G maxn - G[j];880 if (-G[j] >= G maxn2)881 G maxn2 = -G[j];721 double grad_diff = GMaxn - G[j]; 722 if (-G[j] >= GMaxn2) 723 GMaxn2 = -G[j]; 882 724 if (grad_diff > 0) 883 725 { … … 889 731 obj_diff = -(grad_diff * grad_diff) / 1e-12; 890 732 891 if (obj_diff <= obj_diff_ min)733 if (obj_diff <= obj_diff_Min) 892 734 { 893 G min_idx = j;894 obj_diff_ min = obj_diff;735 GMin_idx = j; 736 obj_diff_Min = obj_diff; 895 737 } 896 738 } … … 899 741 } 900 742 901 if (Math.Max(G maxp + Gmaxp2, Gmaxn + Gmaxn2) < eps)743 if (Math.Max(GMaxp + GMaxp2, GMaxn + GMaxn2) < EPS) 902 744 return 1; 903 745 904 if (y[G min_idx] == +1)905 working_set[0] = G maxp_idx;746 if (y[GMin_idx] == +1) 747 working_set[0] = GMaxp_idx; 906 748 else 907 working_set[0] = G maxn_idx;908 working_set[1] = G min_idx;749 working_set[0] = GMaxn_idx; 750 working_set[1] = GMin_idx; 909 751 910 752 return 0; 911 753 } 912 754 913 private bool be_shrunk en(int i, double Gmax1, double Gmax2, double Gmax3, double Gmax4)755 private bool be_shrunk(int i, double GMax1, double GMax2, double GMax3, double GMax4) 914 756 { 915 757 if (is_upper_bound(i)) 916 758 { 917 759 if (y[i] == +1) 918 return (-G[i] > G max1);760 return (-G[i] > GMax1); 919 761 else 920 return (-G[i] > G max4);762 return (-G[i] > GMax4); 921 763 } 922 764 else if (is_lower_bound(i)) 923 765 { 924 766 if (y[i] == +1) 925 return (G[i] > G max2);767 return (G[i] > GMax2); 926 768 else 927 return (G[i] > G max3);769 return (G[i] > GMax3); 928 770 } 929 771 else … … 931 773 } 932 774 933 pr otected override void do_shrinking()934 { 935 double G max1 = -INF; // max { -y_i * grad(f)_i | y_i = +1, i in I_up(\alpha) }936 double G max2 = -INF; // max { y_i * grad(f)_i | y_i = +1, i in I_low(\alpha) }937 double G max3 = -INF; // max { -y_i * grad(f)_i | y_i = -1, i in I_up(\alpha) }938 double G max4 = -INF; // max { y_i * grad(f)_i | y_i = -1, i in I_low(\alpha) }939 940 // find maximal violating pair first775 private void do_shrinking() 776 { 777 double GMax1 = -INF; // Max { -y_i * grad(f)_i | y_i = +1, i in I_up(\alpha) } 778 double GMax2 = -INF; // Max { y_i * grad(f)_i | y_i = +1, i in I_low(\alpha) } 779 double GMax3 = -INF; // Max { -y_i * grad(f)_i | y_i = -1, i in I_up(\alpha) } 780 double GMax4 = -INF; // Max { y_i * grad(f)_i | y_i = -1, i in I_low(\alpha) } 781 782 // find Maximal violating pair first 941 783 int i; 942 784 for (i = 0; i < active_size; i++) … … 946 788 if (y[i] == +1) 947 789 { 948 if (-G[i] > G max1) Gmax1 = -G[i];949 } 950 else if (-G[i] > G max4) Gmax4 = -G[i];790 if (-G[i] > GMax1) GMax1 = -G[i]; 791 } 792 else if (-G[i] > GMax4) GMax4 = -G[i]; 951 793 } 952 794 if (!is_lower_bound(i)) … … 954 796 if (y[i] == +1) 955 797 { 956 if (G[i] > Gmax2) Gmax2 = G[i]; 957 } 958 else if (G[i] > Gmax3) Gmax3 = G[i]; 959 } 960 } 961 962 // shrinking 798 if (G[i] > GMax2) GMax2 = G[i]; 799 } 800 else if (G[i] > GMax3) GMax3 = G[i]; 801 } 802 } 803 804 if (unshrink == false && Math.Max(GMax1 + GMax2, GMax3 + GMax4) <= EPS * 10) 805 { 806 unshrink = true; 807 reconstruct_gradient(); 808 active_size = l; 809 } 963 810 964 811 for (i = 0; i < active_size; i++) 965 if (be_shrunk en(i, Gmax1, Gmax2, Gmax3, Gmax4))812 if (be_shrunk(i, GMax1, GMax2, GMax3, GMax4)) 966 813 { 967 814 active_size--; 968 815 while (active_size > i) 969 816 { 970 if (!be_shrunk en(active_size, Gmax1, Gmax2, Gmax3, Gmax4))817 if (!be_shrunk(active_size, GMax1, GMax2, GMax3, GMax4)) 971 818 { 972 819 swap_index(i, active_size); … … 976 823 } 977 824 } 978 979 if (unshrinked || Math.Max(Gmax1 + Gmax2, Gmax3 + Gmax4) > eps * 10) return; 980 981 unshrinked = true; 982 reconstruct_gradient(); 983 984 for (i = l - 1; i >= active_size; i--) 985 if (!be_shrunken(i, Gmax1, Gmax2, Gmax3, Gmax4)) 986 { 987 while (active_size < i) 988 { 989 if (be_shrunken(active_size, Gmax1, Gmax2, Gmax3, Gmax4)) 990 { 991 swap_index(i, active_size); 992 break; 993 } 994 active_size++; 995 } 996 active_size++; 997 } 998 } 999 1000 protected override double calculate_rho() 825 } 826 827 private double calculate_rho() 1001 828 { 1002 829 int nr_free1 = 0, nr_free2 = 0; … … 1054 881 class SVC_Q : Kernel 1055 882 { 1056 private s hort[] y;883 private sbyte[] y; 1057 884 private Cache cache; 1058 885 private float[] QD; 1059 886 1060 public SVC_Q(Problem prob, Parameter param, s hort[] y_) : base(prob.Count, prob.X, param)1061 { 1062 y = (s hort[])y_.Clone();887 public SVC_Q(Problem prob, Parameter param, sbyte[] y_) : base(prob.Count, prob.X, param) 888 { 889 y = (sbyte[])y_.Clone(); 1063 890 cache = new Cache(prob.Count, (long)(param.CacheSize * (1 << 20))); 1064 891 QD = new float[prob.Count]; 1065 892 for (int i = 0; i < prob.Count; i++) 1066 QD[i] = (float) kernel_function(i, i);1067 } 1068 1069 public override float[] get_Q(int i, int len)1070 { 1071 float[] [] data = new float[1][];1072 int start ;1073 if ((start = cache. get_data(i,data, len)) < len)1074 { 1075 for ( intj = start; j < len; j++)1076 data[ 0][j] = (float)(y[i] * y[j] * kernel_function(i, j));1077 } 1078 return data [0];1079 } 1080 1081 public override float[] get_QD()893 QD[i] = (float)KernelFunction(i, i); 894 } 895 896 public override sealed float[] GetQ(int i, int len) 897 { 898 float[] data = null; 899 int start, j; 900 if ((start = cache.GetData(i, ref data, len)) < len) 901 { 902 for (j = start; j < len; j++) 903 data[j] = (float)(y[i] * y[j] * KernelFunction(i, j)); 904 } 905 return data; 906 } 907 908 public override sealed float[] GetQD() 1082 909 { 1083 910 return QD; 1084 911 } 1085 912 1086 public override void swap_index(int i, int j)1087 { 1088 cache. swap_index(i, j);1089 base. swap_index(i, j);1090 do { short _ = y[i]; y[i] = y[j]; y[j] = _; } while (false);1091 do { float _ = QD[i]; QD[i] = QD[j]; QD[j] = _; } while (false);913 public override sealed void SwapIndex(int i, int j) 914 { 915 cache.SwapIndex(i, j); 916 base.SwapIndex(i, j); 917 y.SwapIndex(i, j); 918 QD.SwapIndex(i, j); 1092 919 } 1093 920 } … … 1098 925 private float[] QD; 1099 926 1100 public ONE_CLASS_Q(Problem prob, Parameter param) : base(prob.Count, prob.X, param)927 public ONE_CLASS_Q(Problem prob, Parameter param) : base(prob.Count, prob.X, param) 1101 928 { 1102 929 cache = new Cache(prob.Count, (long)(param.CacheSize * (1 << 20))); 1103 930 QD = new float[prob.Count]; 1104 931 for (int i = 0; i < prob.Count; i++) 1105 QD[i] = (float) kernel_function(i, i);1106 } 1107 1108 public override float[] get_Q(int i, int len)1109 { 1110 float[] [] data = new float[1][];1111 int start ;1112 if ((start = cache. get_data(i,data, len)) < len)1113 { 1114 for ( intj = start; j < len; j++)1115 data[ 0][j] = (float)kernel_function(i, j);1116 } 1117 return data [0];1118 } 1119 1120 public override float[] get_QD()932 QD[i] = (float)KernelFunction(i, i); 933 } 934 935 public override sealed float[] GetQ(int i, int len) 936 { 937 float[] data = null; 938 int start, j; 939 if ((start = cache.GetData(i, ref data, len)) < len) 940 { 941 for (j = start; j < len; j++) 942 data[j] = (float)KernelFunction(i, j); 943 } 944 return data; 945 } 946 947 public override sealed float[] GetQD() 1121 948 { 1122 949 return QD; 1123 950 } 1124 951 1125 public override void swap_index(int i, int j)1126 { 1127 cache. swap_index(i, j);1128 base. swap_index(i, j);1129 do { float _ = QD[i]; QD[i] = QD[j]; QD[j] = _; } while (false);952 public override sealed void SwapIndex(int i, int j) 953 { 954 cache.SwapIndex(i, j); 955 base.SwapIndex(i, j); 956 QD.SwapIndex(i, j); 1130 957 } 1131 958 } … … 1135 962 private int l; 1136 963 private Cache cache; 1137 private s hort[] sign;964 private sbyte[] sign; 1138 965 private int[] index; 1139 966 private int next_buffer; … … 1141 968 private float[] QD; 1142 969 1143 public SVR_Q(Problem prob, Parameter param) 1144 : base(prob.Count, prob.X, param) 970 public SVR_Q(Problem prob, Parameter param) : base(prob.Count, prob.X, param) 1145 971 { 1146 972 l = prob.Count; 1147 973 cache = new Cache(l, (long)(param.CacheSize * (1 << 20))); 1148 974 QD = new float[2 * l]; 1149 sign = new s hort[2 * l];975 sign = new sbyte[2 * l]; 1150 976 index = new int[2 * l]; 1151 977 for (int k = 0; k < l; k++) … … 1155 981 index[k] = k; 1156 982 index[k + l] = k; 1157 QD[k] = (float) kernel_function(k, k);983 QD[k] = (float)KernelFunction(k, k); 1158 984 QD[k + l] = QD[k]; 1159 985 } … … 1164 990 } 1165 991 1166 public override void swap_index(int i, int j)1167 { 1168 do { short _ = sign[i]; sign[i] = sign[j]; sign[j] = _; } while (false);1169 do { int _ = index[i]; index[i] = index[j]; index[j] = _; } while (false);1170 do { float _ = QD[i]; QD[i] = QD[j]; QD[j] = _; } while (false);1171 } 1172 1173 public override float[] get_Q(int i, int len)1174 { 1175 float[] [] data = new float[1][];1176 int real_i = index[i];1177 if (cache. get_data(real_i,data, l) < l)1178 { 1179 for ( intj = 0; j < l; j++)1180 data[ 0][j] = (float)kernel_function(real_i, j);992 public override sealed void SwapIndex(int i, int j) 993 { 994 sign.SwapIndex(i, j); 995 index.SwapIndex(i, j); 996 QD.SwapIndex(i, j); 997 } 998 999 public override sealed float[] GetQ(int i, int len) 1000 { 1001 float[] data = null; 1002 int j, real_i = index[i]; 1003 if (cache.GetData(real_i, ref data, l) < l) 1004 { 1005 for (j = 0; j < l; j++) 1006 data[j] = (float)KernelFunction(real_i, j); 1181 1007 } 1182 1008 … … 1184 1010 float[] buf = buffer[next_buffer]; 1185 1011 next_buffer = 1 - next_buffer; 1186 s hortsi = sign[i];1187 for ( intj = 0; j < len; j++)1188 buf[j] = si * sign[j] * data[0][index[j]];1012 sbyte si = sign[i]; 1013 for (j = 0; j < len; j++) 1014 buf[j] = (float)si * sign[j] * data[index[j]]; 1189 1015 return buf; 1190 1016 } 1191 1017 1192 public override float[] get_QD()1018 public override sealed float[] GetQD() 1193 1019 { 1194 1020 return QD; … … 1196 1022 } 1197 1023 1198 internal staticclass Procedures1024 internal class Procedures 1199 1025 { 1026 private static bool _verbose; 1027 public static bool IsVerbose 1028 { 1029 get 1030 { 1031 return _verbose; 1032 } 1033 set 1034 { 1035 _verbose = value; 1036 } 1037 } 1200 1038 // 1201 1039 // construct and solve various formulations 1202 1040 // 1041 public const int LIBSVM_VERSION = 289; 1042 1043 public static TextWriter svm_print_string = Console.Out; 1044 1045 public static void info(string s) 1046 { 1047 if(_verbose) 1048 svm_print_string.Write(s); 1049 } 1050 1203 1051 private static void solve_c_svc(Problem prob, Parameter param, 1204 1052 double[] alpha, Solver.SolutionInfo si, … … 1206 1054 { 1207 1055 int l = prob.Count; 1208 double[] minus_ones = new double[l];1209 s hort[] y = new short[l];1056 double[] Minus_ones = new double[l]; 1057 sbyte[] y = new sbyte[l]; 1210 1058 1211 1059 int i; … … 1214 1062 { 1215 1063 alpha[i] = 0; 1216 minus_ones[i] = -1;1064 Minus_ones[i] = -1; 1217 1065 if (prob.Y[i] > 0) y[i] = +1; else y[i] = -1; 1218 1066 } 1219 1067 1220 1068 Solver s = new Solver(); 1221 s.Solve(l, new SVC_Q(prob, param, y), minus_ones, y,1069 s.Solve(l, new SVC_Q(prob, param, y), Minus_ones, y, 1222 1070 alpha, Cp, Cn, param.EPS, si, param.Shrinking); 1223 1071 … … 1227 1075 1228 1076 if (Cp == Cn) 1229 Debug.Write("nu = " + sum_alpha / (Cp * prob.Count) + "\n");1077 Procedures.info("nu = " + sum_alpha / (Cp * prob.Count) + "\n"); 1230 1078 1231 1079 for (i = 0; i < l; i++) … … 1240 1088 double nu = param.Nu; 1241 1089 1242 s hort[] y = new short[l];1090 sbyte[] y = new sbyte[l]; 1243 1091 1244 1092 for (i = 0; i < l; i++) … … 1269 1117 1270 1118 Solver_NU s = new Solver_NU(); 1271 s.Solve(l, new SVC_Q(prob, param, y), zeros, y, 1272 alpha, 1.0, 1.0, param.EPS, si, param.Shrinking); 1119 s.Solve(l, new SVC_Q(prob, param, y), zeros, y, alpha, 1.0, 1.0, param.EPS, si, param.Shrinking); 1273 1120 double r = si.r; 1274 1121 1275 Debug.Write("C = " + 1 / r + "\n");1122 Procedures.info("C = " + 1 / r + "\n"); 1276 1123 1277 1124 for (i = 0; i < l; i++) … … 1285 1132 1286 1133 private static void solve_one_class(Problem prob, Parameter param, 1287 1134 double[] alpha, Solver.SolutionInfo si) 1288 1135 { 1289 1136 int l = prob.Count; 1290 1137 double[] zeros = new double[l]; 1291 s hort[] ones = new short[l];1138 sbyte[] ones = new sbyte[l]; 1292 1139 int i; 1293 1140 … … 1308 1155 1309 1156 Solver s = new Solver(); 1310 s.Solve(l, new ONE_CLASS_Q(prob, param), zeros, ones, 1311 alpha, 1.0, 1.0, param.EPS, si, param.Shrinking); 1312 } 1313 1314 private static void solve_epsilon_svr(Problem prob, Parameter param, 1315 double[] alpha, Solver.SolutionInfo si) 1157 s.Solve(l, new ONE_CLASS_Q(prob, param), zeros, ones, alpha, 1.0, 1.0, param.EPS, si, param.Shrinking); 1158 } 1159 1160 private static void solve_epsilon_svr(Problem prob, Parameter param, double[] alpha, Solver.SolutionInfo si) 1316 1161 { 1317 1162 int l = prob.Count; 1318 1163 double[] alpha2 = new double[2 * l]; 1319 1164 double[] linear_term = new double[2 * l]; 1320 s hort[] y = new short[2 * l];1165 sbyte[] y = new sbyte[2 * l]; 1321 1166 int i; 1322 1167 … … 1333 1178 1334 1179 Solver s = new Solver(); 1335 s.Solve(2 * l, new SVR_Q(prob, param), linear_term, y, 1336 alpha2, param.C, param.C, param.EPS, si, param.Shrinking); 1180 s.Solve(2 * l, new SVR_Q(prob, param), linear_term, y, alpha2, param.C, param.C, param.EPS, si, param.Shrinking); 1337 1181 1338 1182 double sum_alpha = 0; … … 1342 1186 sum_alpha += Math.Abs(alpha[i]); 1343 1187 } 1344 Debug.Write("nu = " + sum_alpha / (param.C * l) + "\n");1188 Procedures.info("nu = " + sum_alpha / (param.C * l) + "\n"); 1345 1189 } 1346 1190 … … 1352 1196 double[] alpha2 = new double[2 * l]; 1353 1197 double[] linear_term = new double[2 * l]; 1354 s hort[] y = new short[2 * l];1198 sbyte[] y = new sbyte[2 * l]; 1355 1199 int i; 1356 1200 … … 1371 1215 s.Solve(2 * l, new SVR_Q(prob, param), linear_term, y, alpha2, C, C, param.EPS, si, param.Shrinking); 1372 1216 1373 Debug.Write("epsilon = " + (-si.r) + "\n");1217 Procedures.info("epsilon = " + (-si.r) + "\n"); 1374 1218 1375 1219 for (i = 0; i < l; i++) … … 1380 1224 // decision_function 1381 1225 // 1382 privateclass decision_function1226 internal class decision_function 1383 1227 { 1384 1228 public double[] alpha; … … 1386 1230 }; 1387 1231 1388 static decision_function svm_train_one( 1389 Problem prob, Parameter param, 1390 double Cp, double Cn) 1232 static decision_function svm_train_one(Problem prob, Parameter param, double Cp, double Cn) 1391 1233 { 1392 1234 double[] alpha = new double[prob.Count]; … … 1411 1253 } 1412 1254 1413 Debug.Write("obj = " + si.obj + ", rho = " + si.rho + "\n");1255 Procedures.info("obj = " + si.obj + ", rho = " + si.rho + "\n"); 1414 1256 1415 1257 // output SVs … … 1435 1277 } 1436 1278 1437 Debug.Write("nSV = " + nSV + ", nBSV = " + nBSV + "\n");1279 Procedures.info("nSV = " + nSV + ", nBSV = " + nBSV + "\n"); 1438 1280 1439 1281 decision_function f = new decision_function(); … … 1455 1297 else prior0 += 1; 1456 1298 1457 int max_iter = 100;// Maximal number of iterations1458 double min_step = 1e-10; // Minimal step taken in line search1459 double sigma = 1e- 3; // For numerically strict PD of Hessian1299 int Max_iter = 100; // Maximal number of iterations 1300 double Min_step = 1e-10; // Minimal step taken in line search 1301 double sigma = 1e-12; // For numerically strict PD of Hessian 1460 1302 double eps = 1e-5; 1461 1303 double hiTarget = (prior1 + 1.0) / (prior1 + 2.0); … … 1480 1322 fval += (t[i] - 1) * fApB + Math.Log(1 + Math.Exp(fApB)); 1481 1323 } 1482 for (iter = 0; iter < max_iter; iter++)1324 for (iter = 0; iter < Max_iter; iter++) 1483 1325 { 1484 1326 // Update Gradient and Hessian (use H' = H + sigma I) … … 1519 1361 1520 1362 1521 stepsize = 1; 1522 while (stepsize >= min_step)1363 stepsize = 1; // Line Search 1364 while (stepsize >= Min_step) 1523 1365 { 1524 1366 newA = A + stepsize * dA; … … 1545 1387 } 1546 1388 1547 if (stepsize < min_step)1548 { 1549 Debug.Write("Line search fails in two-class probability estimates\n");1389 if (stepsize < Min_step) 1390 { 1391 Procedures.info("Line search fails in two-class probability estimates\n"); 1550 1392 break; 1551 1393 } 1552 1394 } 1553 1395 1554 if (iter >= max_iter)1555 Debug.Write("Reaching maximal iterations in two-class probability estimates\n");1396 if (iter >= Max_iter) 1397 Procedures.info("Reaching Maximal iterations in two-class probability estimates\n"); 1556 1398 probAB[0] = A; probAB[1] = B; 1557 1399 } … … 1568 1410 // Method 2 from the multiclass_prob paper by Wu, Lin, and Weng 1569 1411 private static void multiclass_probability(int k, double[,] r, double[] p) 1570 1571 int t,j;1572 int iter = 0, max_iter=Math.Max(100,k);1573 double[,] Q=new double[k,k];1574 double[] Qp= new double[k];1575 double pQp, eps=0.005/k;1576 1577 for (t=0;t<k;t++)1578 1579 p[t]=1.0/k; // Valid if k = 11580 Q[t,t]=0;1581 for (j=0;j<t;j++)1582 1583 Q[t,t]+=r[j,t]*r[j,t];1584 Q[t,j]=Q[j,t];1585 1586 for (j=t+1;j<k;j++)1587 1588 Q[t,t]+=r[j,t]*r[j,t];1589 Q[t,j]=-r[j,t]*r[t,j];1590 1591 1592 for (iter=0;iter<max_iter;iter++)1593 1594 1595 pQp=0;1596 for (t=0;t<k;t++)1597 1598 Qp[t]=0;1599 for (j=0;j<k;j++)1600 Qp[t]+=Q[t,j]*p[j];1601 pQp+=p[t]*Qp[t];1602 1603 double max_error=0;1604 for (t=0;t<k;t++)1605 1606 double error=Math.Abs(Qp[t]-pQp);1607 if (error>max_error)1608 max_error=error;1609 1610 if (max_error<eps) break;1611 1612 for (t=0;t<k;t++)1613 1614 double diff=(-Qp[t]+pQp)/Q[t,t];1615 p[t]+=diff;1616 pQp=(pQp+diff*(diff*Q[t,t]+2*Qp[t]))/(1+diff)/(1+diff);1617 for (j=0;j<k;j++)1618 1619 Qp[j]=(Qp[j]+diff*Q[t,j])/(1+diff);1620 p[j]/=(1+diff);1621 1622 1623 1624 if (iter>=max_iter)1625 Debug.Write("Exceeds max_iter in multiclass_prob\n");1626 1412 { 1413 int t, j; 1414 int iter = 0, Max_iter = Math.Max(100, k); 1415 double[,] Q = new double[k,k]; 1416 double[] Qp = new double[k]; 1417 double pQp, eps = 0.005 / k; 1418 1419 for (t = 0; t < k; t++) 1420 { 1421 p[t] = 1.0 / k; // Valid if k = 1 1422 Q[t,t] = 0; 1423 for (j = 0; j < t; j++) 1424 { 1425 Q[t,t] += r[j,t] * r[j,t]; 1426 Q[t,j] = Q[j,t]; 1427 } 1428 for (j = t + 1; j < k; j++) 1429 { 1430 Q[t,t] += r[j,t] * r[j,t]; 1431 Q[t,j] = -r[j,t] * r[t,j]; 1432 } 1433 } 1434 for (iter = 0; iter < Max_iter; iter++) 1435 { 1436 // stopping condition, recalculate QP,pQP for numerical accuracy 1437 pQp = 0; 1438 for (t = 0; t < k; t++) 1439 { 1440 Qp[t] = 0; 1441 for (j = 0; j < k; j++) 1442 Qp[t] += Q[t,j] * p[j]; 1443 pQp += p[t] * Qp[t]; 1444 } 1445 double Max_error = 0; 1446 for (t = 0; t < k; t++) 1447 { 1448 double error = Math.Abs(Qp[t] - pQp); 1449 if (error > Max_error) 1450 Max_error = error; 1451 } 1452 if (Max_error < eps) break; 1453 1454 for (t = 0; t < k; t++) 1455 { 1456 double diff = (-Qp[t] + pQp) / Q[t,t]; 1457 p[t] += diff; 1458 pQp = (pQp + diff * (diff * Q[t,t] + 2 * Qp[t])) / (1 + diff) / (1 + diff); 1459 for (j = 0; j < k; j++) 1460 { 1461 Qp[j] = (Qp[j] + diff * Q[t,j]) / (1 + diff); 1462 p[j] /= (1 + diff); 1463 } 1464 } 1465 } 1466 if (iter >= Max_iter) 1467 Procedures.info("Exceeds Max_iter in multiclass_prob\n"); 1468 } 1627 1469 1628 1470 // Cross-validation decision values for probability estimates 1629 1471 private static void svm_binary_svc_probability(Problem prob, Parameter param, double Cp, double Cn, double[] probAB) 1630 1472 { 1631 Random rand = new Random();1632 1473 int i; 1633 1474 int nr_fold = 5; … … 1636 1477 1637 1478 // random shuffle 1479 Random rand = new Random(); 1638 1480 for (i = 0; i < prob.Count; i++) perm[i] = i; 1639 1481 for (i = 0; i < prob.Count; i++) … … 1687 1529 subparam.Probability = false; 1688 1530 subparam.C = 1.0; 1689 subparam.WeightCount = 2; 1690 subparam.WeightLabels = new int[2]; 1691 subparam.Weights = new double[2]; 1692 subparam.WeightLabels[0] = +1; 1693 subparam.WeightLabels[1] = -1; 1694 subparam.Weights[0] = Cp; 1695 subparam.Weights[1] = Cn; 1531 subparam.Weights[1] = Cp; 1532 subparam.Weights[-1] = Cn; 1696 1533 Model submodel = svm_train(subprob, subparam); 1697 1534 for (j = begin; j < end; j++) … … 1718 1555 Parameter newparam = (Parameter)param.Clone(); 1719 1556 newparam.Probability = false; 1720 svm_cross_validation(prob, newparam, nr_fold, ymv , null);1557 svm_cross_validation(prob, newparam, nr_fold, ymv); 1721 1558 for (i = 0; i < prob.Count; i++) 1722 1559 { … … 1734 1571 mae += Math.Abs(ymv[i]); 1735 1572 mae /= (prob.Count - count); 1736 Debug.Write("Prob. model for test data: target value = predicted value + z,\nz: Laplace distribution e^(-|z|/sigma)/(2sigma),sigma=" + mae + "\n");1573 Procedures.info("Prob. model for test data: target value = predicted value + z,\nz: Laplace distribution e^(-|z|/sigma)/(2sigma),sigma=" + mae + "\n"); 1737 1574 return mae; 1738 1575 } … … 1743 1580 { 1744 1581 int l = prob.Count; 1745 int max_nr_class = 16;1582 int Max_nr_class = 16; 1746 1583 int nr_class = 0; 1747 int[] label = new int[ max_nr_class];1748 int[] count = new int[ max_nr_class];1584 int[] label = new int[Max_nr_class]; 1585 int[] count = new int[Max_nr_class]; 1749 1586 int[] data_label = new int[l]; 1750 1587 int i; … … 1765 1602 if (j == nr_class) 1766 1603 { 1767 if (nr_class == max_nr_class)1768 { 1769 max_nr_class *= 2;1770 int[] new_data = new int[ max_nr_class];1604 if (nr_class == Max_nr_class) 1605 { 1606 Max_nr_class *= 2; 1607 int[] new_data = new int[Max_nr_class]; 1771 1608 Array.Copy(label, 0, new_data, 0, label.Length); 1772 1609 label = new_data; 1773 new_data = new int[ max_nr_class];1610 new_data = new int[Max_nr_class]; 1774 1611 Array.Copy(count, 0, new_data, 0, count.Length); 1775 1612 count = new_data; … … 1808 1645 model.Parameter = param; 1809 1646 1810 if (param.SvmType == SvmType.ONE_CLASS || 1647 if (param.SvmType == SvmType.ONE_CLASS || 1811 1648 param.SvmType == SvmType.EPSILON_SVR || 1812 1649 param.SvmType == SvmType.NU_SVR) 1813 1650 { 1814 1651 // regression or one-class-svm 1815 model.NumberOfClasses = 2; 1652 model.NumberOfClasses = 2; 1816 1653 model.ClassLabels = null; 1817 1654 model.NumberOfSVPerClass = null; 1818 1655 model.PairwiseProbabilityA = null; model.PairwiseProbabilityB = null; 1819 1656 model.SupportVectorCoefficients = new double[1][]; 1820 1657 1821 1658 if (param.Probability && 1822 1659 (param.SvmType == SvmType.EPSILON_SVR || … … 1834 1671 int i; 1835 1672 for (i = 0; i < prob.Count; i++) 1836 if (Math.Abs(f.alpha[i]) > 0) ++nSV; 1673 if (Math.Abs(f.alpha[i]) > 0) ++nSV; 1837 1674 model.SupportVectorCount = nSV; 1838 1675 model.SupportVectors = new Node[nSV][]; … … 1873 1710 for (i = 0; i < nr_class; i++) 1874 1711 weighted_C[i] = param.C; 1875 for (i = 0; i < param.WeightCount; i++) 1876 { 1877 int j; 1878 for (j = 0; j < nr_class; j++) 1879 if (param.WeightLabels[i] == label[j]) 1880 break; 1881 if (j == nr_class) 1882 Debug.Write("warning: class label " + param.WeightLabels[i] + " specified in weight is not found\n"); 1883 else 1884 weighted_C[j] *= param.Weights[i]; 1712 foreach (int weightedLabel in param.Weights.Keys) 1713 { 1714 int index = Array.IndexOf<int>(label, weightedLabel); 1715 if (index < 0) 1716 Console.Error.WriteLine("warning: class label " + weightedLabel + " specified in weight is not found"); 1717 else weighted_C[index] *= param.Weights[weightedLabel]; 1885 1718 } 1886 1719 … … 1983 1816 } 1984 1817 1985 Debug.Write("Total nSV = " + nnz + "\n");1818 Procedures.info("Total nSV = " + nnz + "\n"); 1986 1819 1987 1820 model.SupportVectorCount = nnz; … … 2029 1862 2030 1863 // Stratified cross validation 2031 public static void svm_cross_validation(Problem prob, Parameter param, int nr_fold, double[] target , Dictionary<int,double>[] confidence)1864 public static void svm_cross_validation(Problem prob, Parameter param, int nr_fold, double[] target) 2032 1865 { 2033 1866 Random rand = new Random(); … … 2131 1964 param.SvmType == SvmType.NU_SVC)) 2132 1965 { 1966 double[] prob_estimates = new double[svm_get_nr_class(submodel)]; 2133 1967 for (j = begin; j < end; j++) 2134 {2135 double[] prob_estimates = new double[svm_get_nr_class(submodel)];2136 1968 target[perm[j]] = svm_predict_probability(submodel, prob.X[perm[j]], prob_estimates); 2137 confidence[perm[j]] = new Dictionary<int, double>();2138 for (int label = 0; label < prob_estimates.Length; label++)2139 confidence[perm[j]][submodel.ClassLabels[label]] = prob_estimates[label];2140 2141 }2142 1969 } 2143 1970 else … … 2171 1998 else 2172 1999 { 2173 Debug.Write("Model doesn't contain information for SVR probability inference\n");2000 Console.Error.WriteLine("Model doesn't contain information for SVR probability inference"); 2174 2001 return 0; 2175 2002 } … … 2185 2012 double sum = 0; 2186 2013 for (int i = 0; i < model.SupportVectorCount; i++) 2187 sum += sv_coef[i] * Kernel. k_function(x, model.SupportVectors[i], model.Parameter);2014 sum += sv_coef[i] * Kernel.KernelFunction(x, model.SupportVectors[i], model.Parameter); 2188 2015 sum -= model.Rho[0]; 2189 2016 dec_values[0] = sum; … … 2197 2024 double[] kvalue = new double[l]; 2198 2025 for (i = 0; i < l; i++) 2199 kvalue[i] = Kernel. k_function(x, model.SupportVectors[i], model.Parameter);2026 kvalue[i] = Kernel.KernelFunction(x, model.SupportVectors[i], model.Parameter); 2200 2027 2201 2028 int[] start = new int[nr_class]; … … 2262 2089 } 2263 2090 2264 int vote_ max_idx = 0;2091 int vote_Max_idx = 0; 2265 2092 for (i = 1; i < nr_class; i++) 2266 if (vote[i] > vote[vote_ max_idx])2267 vote_ max_idx = i;2268 return model.ClassLabels[vote_ max_idx];2093 if (vote[i] > vote[vote_Max_idx]) 2094 vote_Max_idx = i; 2095 return model.ClassLabels[vote_Max_idx]; 2269 2096 } 2270 2097 } 2271 2098 2272 2099 public static double svm_predict_probability(Model model, Node[] x, double[] prob_estimates) 2273 { 2274 if ((model.Parameter.SvmType == SvmType.C_SVC || model.Parameter.SvmType == SvmType.NU_SVC) && 2275 model.PairwiseProbabilityA!=null && model.PairwiseProbabilityB!=null) 2276 { 2277 int i; 2278 int nr_class = model.NumberOfClasses; 2279 double[] dec_values = new double[nr_class*(nr_class-1)/2]; 2280 svm_predict_values(model, x, dec_values); 2281 2282 double min_prob=1e-7; 2283 double[,] pairwise_prob=new double[nr_class,nr_class]; 2284 2285 int k=0; 2286 for(i=0;i<nr_class;i++) 2287 for(int j=i+1;j<nr_class;j++) 2288 { 2289 pairwise_prob[i,j]=Math.Min(Math.Max(sigmoid_predict(dec_values[k],model.PairwiseProbabilityA[k],model.PairwiseProbabilityB[k]),min_prob),1-min_prob); 2290 pairwise_prob[j,i]=1-pairwise_prob[i,j]; 2291 k++; 2292 } 2293 multiclass_probability(nr_class,pairwise_prob,prob_estimates); 2294 2295 int prob_max_idx = 0; 2296 for(i=1;i<nr_class;i++) 2297 if(prob_estimates[i] > prob_estimates[prob_max_idx]) 2298 prob_max_idx = i; 2299 return model.ClassLabels[prob_max_idx]; 2300 } 2301 else 2302 return svm_predict(model, x); 2303 } 2304 2305 private static double atof(string s) 2306 { 2307 return double.Parse(s); 2308 } 2309 2310 private static int atoi(string s) 2311 { 2312 return int.Parse(s); 2100 { 2101 if ((model.Parameter.SvmType == SvmType.C_SVC || model.Parameter.SvmType == SvmType.NU_SVC) && 2102 model.PairwiseProbabilityA != null && model.PairwiseProbabilityB != null) 2103 { 2104 int i; 2105 int nr_class = model.NumberOfClasses; 2106 double[] dec_values = new double[nr_class * (nr_class - 1) / 2]; 2107 svm_predict_values(model, x, dec_values); 2108 2109 double Min_prob = 1e-7; 2110 double[,] pairwise_prob = new double[nr_class, nr_class]; 2111 2112 int k = 0; 2113 for (i = 0; i < nr_class; i++) 2114 { 2115 for (int j = i + 1; j < nr_class; j++) 2116 { 2117 pairwise_prob[i, j] = Math.Min(Math.Max(sigmoid_predict(dec_values[k], model.PairwiseProbabilityA[k], model.PairwiseProbabilityB[k]), Min_prob), 1 - Min_prob); 2118 pairwise_prob[j, i] = 1 - pairwise_prob[i, j]; 2119 k++; 2120 } 2121 } 2122 multiclass_probability(nr_class, pairwise_prob, prob_estimates); 2123 2124 int prob_Max_idx = 0; 2125 for (i = 1; i < nr_class; i++) 2126 if (prob_estimates[i] > prob_estimates[prob_Max_idx]) 2127 prob_Max_idx = i; 2128 return model.ClassLabels[prob_Max_idx]; 2129 } 2130 else 2131 return svm_predict(model, x); 2313 2132 } 2314 2133 … … 2318 2137 2319 2138 SvmType svm_type = param.SvmType; 2320 if (svm_type != SvmType.C_SVC &&2321 svm_type != SvmType.NU_SVC &&2322 svm_type != SvmType.ONE_CLASS &&2323 svm_type != SvmType.EPSILON_SVR &&2324 svm_type != SvmType.NU_SVR)2325 return "unknown svm type";2326 2139 2327 2140 // kernel_type, degree 2328 2141 2329 2142 KernelType kernel_type = param.KernelType; 2330 if (kernel_type != KernelType.LINEAR &&2331 kernel_type != KernelType.POLY &&2332 kernel_type != KernelType.RBF &&2333 kernel_type != KernelType.SIGMOID &&2334 kernel_type != KernelType.PRECOMPUTED)2335 return "unknown kernel type";2336 2143 2337 2144 if (param.Degree < 0) … … 2345 2152 if (param.EPS <= 0) 2346 2153 return "eps <= 0"; 2154 2155 if (param.Gamma == 0) 2156 param.Gamma = 1.0 / prob.MaxIndex; 2347 2157 2348 2158 if (svm_type == SvmType.C_SVC || … … 2362 2172 return "p < 0"; 2363 2173 2364 if (param.Probability && svm_type == SvmType.ONE_CLASS) 2174 if (param.Probability && 2175 svm_type == SvmType.ONE_CLASS) 2365 2176 return "one-class SVM probability output not supported yet"; 2366 2177 … … 2370 2181 { 2371 2182 int l = prob.Count; 2372 int max_nr_class = 16;2183 int Max_nr_class = 16; 2373 2184 int nr_class = 0; 2374 int[] label = new int[ max_nr_class];2375 int[] count = new int[ max_nr_class];2185 int[] label = new int[Max_nr_class]; 2186 int[] count = new int[Max_nr_class]; 2376 2187 2377 2188 int i; … … 2389 2200 if (j == nr_class) 2390 2201 { 2391 if (nr_class == max_nr_class)2392 { 2393 max_nr_class *= 2;2394 int[] new_data = new int[ max_nr_class];2202 if (nr_class == Max_nr_class) 2203 { 2204 Max_nr_class *= 2; 2205 int[] new_data = new int[Max_nr_class]; 2395 2206 Array.Copy(label, 0, new_data, 0, label.Length); 2396 2207 label = new_data; 2397 2208 2398 new_data = new int[ max_nr_class];2209 new_data = new int[Max_nr_class]; 2399 2210 Array.Copy(count, 0, new_data, 0, count.Length); 2400 2211 count = new_data; … … 2432 2243 } 2433 2244 } 2434 2435 2245 }
Note: See TracChangeset
for help on using the changeset viewer.