1  /*


2  * SVM.NET Library


3  * Copyright (C) 2008 Matthew Johnson


4  *


5  * This program is free software: you can redistribute it and/or modify


6  * it under the terms of the GNU General Public License as published by


7  * the Free Software Foundation, either version 3 of the License, or


8  * (at your option) any later version.


9  *


10  * This program is distributed in the hope that it will be useful,


11  * but WITHOUT ANY WARRANTY; without even the implied warranty of


12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the


13  * GNU General Public License for more details.


14  *


15  * You should have received a copy of the GNU General Public License


16  * along with this program. If not, see <http://www.gnu.org/licenses/>.


17  */


18 


19 


20  using System;


21  using System.IO;


22 


23  namespace SVM {


24 


25  // An SMO algorithm in Fan et al., JMLR 6(2005), p. 18891918


26  // Solves:


27  //


28  // Min 0.5(\alpha^T Q \alpha) + p^T \alpha


29  //


30  // y^T \alpha = \delta


31  // y_i = +1 or 1


32  // 0 <= alpha_i <= Cp for y_i = 1


33  // 0 <= alpha_i <= Cn for y_i = 1


34  //


35  // Given:


36  //


37  // Q, p, y, Cp, Cn, and an initial feasible point \alpha


38  // l is the size of vectors and matrices


39  // eps is the stopping tolerance


40  //


41  // solution will be put in \alpha, objective value will be put in obj


42  //


43  internal class Solver {


44  protected int active_size;


45  protected sbyte[] y;


46  protected double[] G; // gradient of objective function


47  private const byte LOWER_BOUND = 0;


48  private const byte UPPER_BOUND = 1;


49  private const byte FREE = 2;


50  private byte[] alpha_status; // LOWER_BOUND, UPPER_BOUND, FREE


51  private double[] alpha;


52  protected IQMatrix Q;


53  protected float[] QD;


54  protected double EPS;


55  private double Cp, Cn;


56  private double[] p;


57  private int[] active_set;


58  private double[] G_bar; // gradient, if we treat free variables as 0


59  protected int l;


60  protected bool unshrink; // XXX


61 


62  protected const double INF = double.PositiveInfinity;


63 


64  private double get_C(int i) {


65  return (y[i] > 0) ? Cp : Cn;


66  }


67 


68  private void update_alpha_status(int i) {


69  if (alpha[i] >= get_C(i))


70  alpha_status[i] = UPPER_BOUND;


71  else if (alpha[i] <= 0)


72  alpha_status[i] = LOWER_BOUND;


73  else alpha_status[i] = FREE;


74  }


75 


76  protected bool is_upper_bound(int i) { return alpha_status[i] == UPPER_BOUND; }


77  protected bool is_lower_bound(int i) { return alpha_status[i] == LOWER_BOUND; }


78 


79  private bool is_free(int i) { return alpha_status[i] == FREE; }


80 


81  public class SolutionInfo {


82  public double obj;


83  public double rho;


84  public double upper_bound_p;


85  public double upper_bound_n;


86  public double r; // for Solver_NU


87  }


88 


89  protected void swap_index(int i, int j) {


90  Q.SwapIndex(i, j);


91  y.SwapIndex(i, j);


92  G.SwapIndex(i, j);


93  alpha_status.SwapIndex(i, j);


94  alpha.SwapIndex(i, j);


95  p.SwapIndex(i, j);


96  active_set.SwapIndex(i, j);


97  G_bar.SwapIndex(i, j);


98  }


99 


100  protected void reconstruct_gradient() {


101  // reconstruct inactive elements of G from G_bar and free variables


102 


103  if (active_size == l) return;


104 


105  int i, j;


106  int nr_free = 0;


107 


108  for (j = active_size; j < l; j++)


109  G[j] = G_bar[j] + p[j];


110 


111  for (j = 0; j < active_size; j++)


112  if (is_free(j))


113  nr_free++;


114 


115  if (2 * nr_free < active_size)


116  Procedures.info("\nWarning: using h 0 may be faster\n");


117 


118  if (nr_free * l > 2 * active_size * (l  active_size)) {


119  for (i = active_size; i < l; i++) {


120  float[] Q_i = Q.GetQ(i, active_size);


121  for (j = 0; j < active_size; j++)


122  if (is_free(j))


123  G[i] += alpha[j] * Q_i[j];


124  }


125  } else {


126  for (i = 0; i < active_size; i++)


127  if (is_free(i)) {


128  float[] Q_i = Q.GetQ(i, l);


129  double alpha_i = alpha[i];


130  for (j = active_size; j < l; j++)


131  G[j] += alpha_i * Q_i[j];


132  }


133  }


134  }


135 


136  public virtual void Solve(int l, IQMatrix Q, double[] p_, sbyte[] y_, double[] alpha_, double Cp, double Cn, double eps, SolutionInfo si, bool shrinking) {


137  this.l = l;


138  this.Q = Q;


139  QD = Q.GetQD();


140  p = (double[])p_.Clone();


141  y = (sbyte[])y_.Clone();


142  alpha = (double[])alpha_.Clone();


143  this.Cp = Cp;


144  this.Cn = Cn;


145  this.EPS = eps;


146  this.unshrink = false;


147 


148  // initialize alpha_status


149  {


150  alpha_status = new byte[l];


151  for (int i = 0; i < l; i++)


152  update_alpha_status(i);


153  }


154 


155  // initialize active set (for shrinking)


156  {


157  active_set = new int[l];


158  for (int i = 0; i < l; i++)


159  active_set[i] = i;


160  active_size = l;


161  }


162 


163  // initialize gradient


164  {


165  G = new double[l];


166  G_bar = new double[l];


167  int i;


168  for (i = 0; i < l; i++) {


169  G[i] = p[i];


170  G_bar[i] = 0;


171  }


172  for (i = 0; i < l; i++)


173  if (!is_lower_bound(i)) {


174  float[] Q_i = Q.GetQ(i, l);


175  double alpha_i = alpha[i];


176  int j;


177  for (j = 0; j < l; j++)


178  G[j] += alpha_i * Q_i[j];


179  if (is_upper_bound(i))


180  for (j = 0; j < l; j++)


181  G_bar[j] += get_C(i) * Q_i[j];


182  }


183  }


184 


185  // optimization step


186 


187  int iter = 0;


188  int counter = Math.Min(l, 1000) + 1;


189  int[] working_set = new int[2];


190 


191  while (true) {


192  // show progress and do shrinking


193 


194  if (counter == 0) {


195  counter = Math.Min(l, 1000);


196  if (shrinking) do_shrinking();


197  Procedures.info(".");


198  }


199 


200  if (select_working_set(working_set) != 0) {


201  // reconstruct the whole gradient


202  reconstruct_gradient();


203  // reset active set size and check


204  active_size = l;


205  Procedures.info("*");


206  if (select_working_set(working_set) != 0)


207  break;


208  else


209  counter = 1; // do shrinking next iteration


210  }


211 


212  int i = working_set[0];


213  int j = working_set[1];


214 


215  ++iter;


216 


217  // update alpha[i] and alpha[j], handle bounds carefully


218 


219  float[] Q_i = Q.GetQ(i, active_size);


220  float[] Q_j = Q.GetQ(j, active_size);


221 


222  double C_i = get_C(i);


223  double C_j = get_C(j);


224 


225  double old_alpha_i = alpha[i];


226  double old_alpha_j = alpha[j];


227 


228  if (y[i] != y[j]) {


229  double quad_coef = Q_i[i] + Q_j[j] + 2 * Q_i[j];


230  if (quad_coef <= 0)


231  quad_coef = 1e12;


232  double delta = (G[i]  G[j]) / quad_coef;


233  double diff = alpha[i]  alpha[j];


234  alpha[i] += delta;


235  alpha[j] += delta;


236 


237  if (diff > 0) {


238  if (alpha[j] < 0) {


239  alpha[j] = 0;


240  alpha[i] = diff;


241  }


242  } else {


243  if (alpha[i] < 0) {


244  alpha[i] = 0;


245  alpha[j] = diff;


246  }


247  }


248  if (diff > C_i  C_j) {


249  if (alpha[i] > C_i) {


250  alpha[i] = C_i;


251  alpha[j] = C_i  diff;


252  }


253  } else {


254  if (alpha[j] > C_j) {


255  alpha[j] = C_j;


256  alpha[i] = C_j + diff;


257  }


258  }


259  } else {


260  double quad_coef = Q_i[i] + Q_j[j]  2 * Q_i[j];


261  if (quad_coef <= 0)


262  quad_coef = 1e12;


263  double delta = (G[i]  G[j]) / quad_coef;


264  double sum = alpha[i] + alpha[j];


265  alpha[i] = delta;


266  alpha[j] += delta;


267 


268  if (sum > C_i) {


269  if (alpha[i] > C_i) {


270  alpha[i] = C_i;


271  alpha[j] = sum  C_i;


272  }


273  } else {


274  if (alpha[j] < 0) {


275  alpha[j] = 0;


276  alpha[i] = sum;


277  }


278  }


279  if (sum > C_j) {


280  if (alpha[j] > C_j) {


281  alpha[j] = C_j;


282  alpha[i] = sum  C_j;


283  }


284  } else {


285  if (alpha[i] < 0) {


286  alpha[i] = 0;


287  alpha[j] = sum;


288  }


289  }


290  }


291 


292  // update G


293 


294  double delta_alpha_i = alpha[i]  old_alpha_i;


295  double delta_alpha_j = alpha[j]  old_alpha_j;


296 


297  for (int k = 0; k < active_size; k++) {


298  G[k] += Q_i[k] * delta_alpha_i + Q_j[k] * delta_alpha_j;


299  }


300 


301  // update alpha_status and G_bar


302 


303  {


304  bool ui = is_upper_bound(i);


305  bool uj = is_upper_bound(j);


306  update_alpha_status(i);


307  update_alpha_status(j);


308  int k;


309  if (ui != is_upper_bound(i)) {


310  Q_i = Q.GetQ(i, l);


311  if (ui)


312  for (k = 0; k < l; k++)


313  G_bar[k] = C_i * Q_i[k];


314  else


315  for (k = 0; k < l; k++)


316  G_bar[k] += C_i * Q_i[k];


317  }


318 


319  if (uj != is_upper_bound(j)) {


320  Q_j = Q.GetQ(j, l);


321  if (uj)


322  for (k = 0; k < l; k++)


323  G_bar[k] = C_j * Q_j[k];


324  else


325  for (k = 0; k < l; k++)


326  G_bar[k] += C_j * Q_j[k];


327  }


328  }


329 


330  }


331 


332  // calculate rho


333 


334  si.rho = calculate_rho();


335 


336  // calculate objective value


337  {


338  double v = 0;


339  int i;


340  for (i = 0; i < l; i++)


341  v += alpha[i] * (G[i] + p[i]);


342 


343  si.obj = v / 2;


344  }


345 


346  // put back the solution


347  {


348  for (int i = 0; i < l; i++)


349  alpha_[active_set[i]] = alpha[i];


350  }


351 


352  si.upper_bound_p = Cp;


353  si.upper_bound_n = Cn;


354 


355  Procedures.info("\noptimization finished, #iter = " + iter + "\n");


356  }


357 


358  // return 1 if already optimal, return 0 otherwise


359  protected virtual int select_working_set(int[] working_set) {


360  // return i,j such that


361  // i: Maximizes y_i * grad(f)_i, i in I_up(\alpha)


362  // j: mimimizes the decrease of obj value


363  // (if quadratic coefficeint <= 0, replace it with tau)


364  // y_j*grad(f)_j < y_i*grad(f)_i, j in I_low(\alpha)


365 


366  double GMax = INF;


367  double GMax2 = INF;


368  int GMax_idx = 1;


369  int GMin_idx = 1;


370  double obj_diff_Min = INF;


371 


372  for (int t = 0; t < active_size; t++)


373  if (y[t] == +1) {


374  if (!is_upper_bound(t))


375  if (G[t] >= GMax) {


376  GMax = G[t];


377  GMax_idx = t;


378  }


379  } else {


380  if (!is_lower_bound(t))


381  if (G[t] >= GMax) {


382  GMax = G[t];


383  GMax_idx = t;


384  }


385  }


386 


387  int i = GMax_idx;


388  float[] Q_i = null;


389  if (i != 1) // null Q_i not accessed: GMax=INF if i=1


390  Q_i = Q.GetQ(i, active_size);


391 


392  for (int j = 0; j < active_size; j++) {


393  if (y[j] == +1) {


394  if (!is_lower_bound(j)) {


395  double grad_diff = GMax + G[j];


396  if (G[j] >= GMax2)


397  GMax2 = G[j];


398  if (grad_diff > 0) {


399  double obj_diff;


400  double quad_coef = Q_i[i] + QD[j]  2.0 * y[i] * Q_i[j];


401  if (quad_coef > 0)


402  obj_diff = (grad_diff * grad_diff) / quad_coef;


403  else


404  obj_diff = (grad_diff * grad_diff) / 1e12;


405 


406  if (obj_diff <= obj_diff_Min) {


407  GMin_idx = j;


408  obj_diff_Min = obj_diff;


409  }


410  }


411  }


412  } else {


413  if (!is_upper_bound(j)) {


414  double grad_diff = GMax  G[j];


415  if (G[j] >= GMax2)


416  GMax2 = G[j];


417  if (grad_diff > 0) {


418  double obj_diff;


419  double quad_coef = Q_i[i] + QD[j] + 2.0 * y[i] * Q_i[j];


420  if (quad_coef > 0)


421  obj_diff = (grad_diff * grad_diff) / quad_coef;


422  else


423  obj_diff = (grad_diff * grad_diff) / 1e12;


424 


425  if (obj_diff <= obj_diff_Min) {


426  GMin_idx = j;


427  obj_diff_Min = obj_diff;


428  }


429  }


430  }


431  }


432  }


433 


434  if (GMax + GMax2 < EPS)


435  return 1;


436 


437  working_set[0] = GMax_idx;


438  working_set[1] = GMin_idx;


439  return 0;


440  }


441 


442  private bool be_shrunk(int i, double GMax1, double GMax2) {


443  if (is_upper_bound(i)) {


444  if (y[i] == +1)


445  return (G[i] > GMax1);


446  else


447  return (G[i] > GMax2);


448  } else if (is_lower_bound(i)) {


449  if (y[i] == +1)


450  return (G[i] > GMax2);


451  else


452  return (G[i] > GMax1);


453  } else


454  return (false);


455  }


456 


457  protected virtual void do_shrinking() {


458  int i;


459  double GMax1 = INF; // Max { y_i * grad(f)_i  i in I_up(\alpha) }


460  double GMax2 = INF; // Max { y_i * grad(f)_i  i in I_low(\alpha) }


461 


462  // find Maximal violating pair first


463  for (i = 0; i < active_size; i++) {


464  if (y[i] == +1) {


465  if (!is_upper_bound(i)) {


466  if (G[i] >= GMax1)


467  GMax1 = G[i];


468  }


469  if (!is_lower_bound(i)) {


470  if (G[i] >= GMax2)


471  GMax2 = G[i];


472  }


473  } else {


474  if (!is_upper_bound(i)) {


475  if (G[i] >= GMax2)


476  GMax2 = G[i];


477  }


478  if (!is_lower_bound(i)) {


479  if (G[i] >= GMax1)


480  GMax1 = G[i];


481  }


482  }


483  }


484 


485  if (unshrink == false && GMax1 + GMax2 <= EPS * 10) {


486  unshrink = true;


487  reconstruct_gradient();


488  active_size = l;


489  }


490 


491  for (i = 0; i < active_size; i++)


492  if (be_shrunk(i, GMax1, GMax2)) {


493  active_size;


494  while (active_size > i) {


495  if (!be_shrunk(active_size, GMax1, GMax2)) {


496  swap_index(i, active_size);


497  break;


498  }


499  active_size;


500  }


501  }


502  }


503 


504  protected virtual double calculate_rho() {


505  double r;


506  int nr_free = 0;


507  double ub = INF, lb = INF, sum_free = 0;


508  for (int i = 0; i < active_size; i++) {


509  double yG = y[i] * G[i];


510 


511  if (is_lower_bound(i)) {


512  if (y[i] > 0)


513  ub = Math.Min(ub, yG);


514  else


515  lb = Math.Max(lb, yG);


516  } else if (is_upper_bound(i)) {


517  if (y[i] < 0)


518  ub = Math.Min(ub, yG);


519  else


520  lb = Math.Max(lb, yG);


521  } else {


522  ++nr_free;


523  sum_free += yG;


524  }


525  }


526 


527  if (nr_free > 0)


528  r = sum_free / nr_free;


529  else


530  r = (ub + lb) / 2;


531 


532  return r;


533  }


534 


535  }


536 


537  //


538  // Solver for nusvm classification and regression


539  //


540  // additional constraint: e^T \alpha = constant


541  //


542  class Solver_NU : Solver {


543  private SolutionInfo si;


544 


545  public sealed override void Solve(int l, IQMatrix Q, double[] p, sbyte[] y,


546  double[] alpha, double Cp, double Cn, double eps,


547  SolutionInfo si, bool shrinking) {


548  this.si = si;


549  base.Solve(l, Q, p, y, alpha, Cp, Cn, eps, si, shrinking);


550  }


551 


552  // return 1 if already optimal, return 0 otherwise


553  protected override sealed int select_working_set(int[] working_set) {


554  // return i,j such that y_i = y_j and


555  // i: Maximizes y_i * grad(f)_i, i in I_up(\alpha)


556  // j: Minimizes the decrease of obj value


557  // (if quadratic coefficeint <= 0, replace it with tau)


558  // y_j*grad(f)_j < y_i*grad(f)_i, j in I_low(\alpha)


559 


560  double GMaxp = INF;


561  double GMaxp2 = INF;


562  int GMaxp_idx = 1;


563 


564  double GMaxn = INF;


565  double GMaxn2 = INF;


566  int GMaxn_idx = 1;


567 


568  int GMin_idx = 1;


569  double obj_diff_Min = INF;


570 


571  for (int t = 0; t < active_size; t++)


572  if (y[t] == +1) {


573  if (!is_upper_bound(t))


574  if (G[t] >= GMaxp) {


575  GMaxp = G[t];


576  GMaxp_idx = t;


577  }


578  } else {


579  if (!is_lower_bound(t))


580  if (G[t] >= GMaxn) {


581  GMaxn = G[t];


582  GMaxn_idx = t;


583  }


584  }


585 


586  int ip = GMaxp_idx;


587  int iN = GMaxn_idx;


588  float[] Q_ip = null;


589  float[] Q_in = null;


590  if (ip != 1) // null Q_ip not accessed: GMaxp=INF if ip=1


591  Q_ip = Q.GetQ(ip, active_size);


592  if (iN != 1)


593  Q_in = Q.GetQ(iN, active_size);


594 


595  for (int j = 0; j < active_size; j++) {


596  if (y[j] == +1) {


597  if (!is_lower_bound(j)) {


598  double grad_diff = GMaxp + G[j];


599  if (G[j] >= GMaxp2)


600  GMaxp2 = G[j];


601  if (grad_diff > 0) {


602  double obj_diff;


603  double quad_coef = Q_ip[ip] + QD[j]  2 * Q_ip[j];


604  if (quad_coef > 0)


605  obj_diff = (grad_diff * grad_diff) / quad_coef;


606  else


607  obj_diff = (grad_diff * grad_diff) / 1e12;


608 


609  if (obj_diff <= obj_diff_Min) {


610  GMin_idx = j;


611  obj_diff_Min = obj_diff;


612  }


613  }


614  }


615  } else {


616  if (!is_upper_bound(j)) {


617  double grad_diff = GMaxn  G[j];


618  if (G[j] >= GMaxn2)


619  GMaxn2 = G[j];


620  if (grad_diff > 0) {


621  double obj_diff;


622  double quad_coef = Q_in[iN] + QD[j]  2 * Q_in[j];


623  if (quad_coef > 0)


624  obj_diff = (grad_diff * grad_diff) / quad_coef;


625  else


626  obj_diff = (grad_diff * grad_diff) / 1e12;


627 


628  if (obj_diff <= obj_diff_Min) {


629  GMin_idx = j;


630  obj_diff_Min = obj_diff;


631  }


632  }


633  }


634  }


635  }


636 


637  if (Math.Max(GMaxp + GMaxp2, GMaxn + GMaxn2) < EPS)


638  return 1;


639 


640  if (y[GMin_idx] == +1)


641  working_set[0] = GMaxp_idx;


642  else


643  working_set[0] = GMaxn_idx;


644  working_set[1] = GMin_idx;


645 


646  return 0;


647  }


648 


649  private bool be_shrunk(int i, double GMax1, double GMax2, double GMax3, double GMax4) {


650  if (is_upper_bound(i)) {


651  if (y[i] == +1)


652  return (G[i] > GMax1);


653  else


654  return (G[i] > GMax4);


655  } else if (is_lower_bound(i)) {


656  if (y[i] == +1)


657  return (G[i] > GMax2);


658  else


659  return (G[i] > GMax3);


660  } else


661  return (false);


662  }


663 


664  protected override sealed void do_shrinking() {


665  double GMax1 = INF; // Max { y_i * grad(f)_i  y_i = +1, i in I_up(\alpha) }


666  double GMax2 = INF; // Max { y_i * grad(f)_i  y_i = +1, i in I_low(\alpha) }


667  double GMax3 = INF; // Max { y_i * grad(f)_i  y_i = 1, i in I_up(\alpha) }


668  double GMax4 = INF; // Max { y_i * grad(f)_i  y_i = 1, i in I_low(\alpha) }


669 


670  // find Maximal violating pair first


671  int i;


672  for (i = 0; i < active_size; i++) {


673  if (!is_upper_bound(i)) {


674  if (y[i] == +1) {


675  if (G[i] > GMax1) GMax1 = G[i];


676  } else if (G[i] > GMax4) GMax4 = G[i];


677  }


678  if (!is_lower_bound(i)) {


679  if (y[i] == +1) {


680  if (G[i] > GMax2) GMax2 = G[i];


681  } else if (G[i] > GMax3) GMax3 = G[i];


682  }


683  }


684 


685  if (unshrink == false && Math.Max(GMax1 + GMax2, GMax3 + GMax4) <= EPS * 10) {


686  unshrink = true;


687  reconstruct_gradient();


688  active_size = l;


689  }


690 


691  for (i = 0; i < active_size; i++)


692  if (be_shrunk(i, GMax1, GMax2, GMax3, GMax4)) {


693  active_size;


694  while (active_size > i) {


695  if (!be_shrunk(active_size, GMax1, GMax2, GMax3, GMax4)) {


696  swap_index(i, active_size);


697  break;


698  }


699  active_size;


700  }


701  }


702  }


703 


704  protected override sealed double calculate_rho() {


705  int nr_free1 = 0, nr_free2 = 0;


706  double ub1 = INF, ub2 = INF;


707  double lb1 = INF, lb2 = INF;


708  double sum_free1 = 0, sum_free2 = 0;


709 


710  for (int i = 0; i < active_size; i++) {


711  if (y[i] == +1) {


712  if (is_lower_bound(i))


713  ub1 = Math.Min(ub1, G[i]);


714  else if (is_upper_bound(i))


715  lb1 = Math.Max(lb1, G[i]);


716  else {


717  ++nr_free1;


718  sum_free1 += G[i];


719  }


720  } else {


721  if (is_lower_bound(i))


722  ub2 = Math.Min(ub2, G[i]);


723  else if (is_upper_bound(i))


724  lb2 = Math.Max(lb2, G[i]);


725  else {


726  ++nr_free2;


727  sum_free2 += G[i];


728  }


729  }


730  }


731 


732  double r1, r2;


733  if (nr_free1 > 0)


734  r1 = sum_free1 / nr_free1;


735  else


736  r1 = (ub1 + lb1) / 2;


737 


738  if (nr_free2 > 0)


739  r2 = sum_free2 / nr_free2;


740  else


741  r2 = (ub2 + lb2) / 2;


742 


743  si.r = (r1 + r2) / 2;


744  return (r1  r2) / 2;


745  }


746  }


747 


748  //


749  // Q matrices for various formulations


750  //


751  class SVC_Q : Kernel {


752  private sbyte[] y;


753  private Cache cache;


754  private float[] QD;


755 


756  public SVC_Q(Problem prob, Parameter param, sbyte[] y_)


757  : base(prob.Count, prob.X, param) {


758  y = (sbyte[])y_.Clone();


759  cache = new Cache(prob.Count, (long)(param.CacheSize * (1 << 20)));


760  QD = new float[prob.Count];


761  for (int i = 0; i < prob.Count; i++)


762  QD[i] = (float)KernelFunction(i, i);


763  }


764 


765  public override sealed float[] GetQ(int i, int len) {


766  float[] data = null;


767  int start, j;


768  if ((start = cache.GetData(i, ref data, len)) < len) {


769  for (j = start; j < len; j++)


770  data[j] = (float)(y[i] * y[j] * KernelFunction(i, j));


771  }


772  return data;


773  }


774 


775  public override sealed float[] GetQD() {


776  return QD;


777  }


778 


779  public override sealed void SwapIndex(int i, int j) {


780  cache.SwapIndex(i, j);


781  base.SwapIndex(i, j);


782  y.SwapIndex(i, j);


783  QD.SwapIndex(i, j);


784  }


785  }


786 


787  class ONE_CLASS_Q : Kernel {


788  private Cache cache;


789  private float[] QD;


790 


791  public ONE_CLASS_Q(Problem prob, Parameter param)


792  : base(prob.Count, prob.X, param) {


793  cache = new Cache(prob.Count, (long)(param.CacheSize * (1 << 20)));


794  QD = new float[prob.Count];


795  for (int i = 0; i < prob.Count; i++)


796  QD[i] = (float)KernelFunction(i, i);


797  }


798 


799  public override sealed float[] GetQ(int i, int len) {


800  float[] data = null;


801  int start, j;


802  if ((start = cache.GetData(i, ref data, len)) < len) {


803  for (j = start; j < len; j++)


804  data[j] = (float)KernelFunction(i, j);


805  }


806  return data;


807  }


808 


809  public override sealed float[] GetQD() {


810  return QD;


811  }


812 


813  public override sealed void SwapIndex(int i, int j) {


814  cache.SwapIndex(i, j);


815  base.SwapIndex(i, j);


816  QD.SwapIndex(i, j);


817  }


818  }


819 


820  class SVR_Q : Kernel {


821  private int l;


822  private Cache cache;


823  private sbyte[] sign;


824  private int[] index;


825  private int next_buffer;


826  private float[][] buffer;


827  private float[] QD;


828 


829  public SVR_Q(Problem prob, Parameter param)


830  : base(prob.Count, prob.X, param) {


831  l = prob.Count;


832  cache = new Cache(l, (long)(param.CacheSize * (1 << 20)));


833  QD = new float[2 * l];


834  sign = new sbyte[2 * l];


835  index = new int[2 * l];


836  for (int k = 0; k < l; k++) {


837  sign[k] = 1;


838  sign[k + l] = 1;


839  index[k] = k;


840  index[k + l] = k;


841  QD[k] = (float)KernelFunction(k, k);


842  QD[k + l] = QD[k];


843  }


844  buffer = new float[2][];


845  buffer[0] = new float[2 * l];


846  buffer[1] = new float[2 * l];


847  next_buffer = 0;


848  }


849 


850  public override sealed void SwapIndex(int i, int j) {


851  sign.SwapIndex(i, j);


852  index.SwapIndex(i, j);


853  QD.SwapIndex(i, j);


854  }


855 


856  public override sealed float[] GetQ(int i, int len) {


857  float[] data = null;


858  int j, real_i = index[i];


859  if (cache.GetData(real_i, ref data, l) < l) {


860  for (j = 0; j < l; j++)


861  data[j] = (float)KernelFunction(real_i, j);


862  }


863 


864  // reorder and copy


865  float[] buf = buffer[next_buffer];


866  next_buffer = 1  next_buffer;


867  sbyte si = sign[i];


868  for (j = 0; j < len; j++)


869  buf[j] = (float)si * sign[j] * data[index[j]];


870  return buf;


871  }


872 


873  public override sealed float[] GetQD() {


874  return QD;


875  }


876  }


877 


878  internal class Procedures {


879  private static bool _verbose;


880  public static bool IsVerbose {


881  get {


882  return _verbose;


883  }


884  set {


885  _verbose = value;


886  }


887  }


888  //


889  // construct and solve various formulations


890  //


891  public const int LIBSVM_VERSION = 289;


892 


893  public static TextWriter svm_print_string = Console.Out;


894 


895  public static void info(string s) {


896  if (_verbose)


897  svm_print_string.Write(s);


898  }


899 


900  private static void solve_c_svc(Problem prob, Parameter param,


901  double[] alpha, Solver.SolutionInfo si,


902  double Cp, double Cn) {


903  int l = prob.Count;


904  double[] Minus_ones = new double[l];


905  sbyte[] y = new sbyte[l];


906 


907  int i;


908 


909  for (i = 0; i < l; i++) {


910  alpha[i] = 0;


911  Minus_ones[i] = 1;


912  if (prob.Y[i] > 0) y[i] = +1; else y[i] = 1;


913  }


914 


915  Solver s = new Solver();


916  s.Solve(l, new SVC_Q(prob, param, y), Minus_ones, y,


917  alpha, Cp, Cn, param.EPS, si, param.Shrinking);


918 


919  double sum_alpha = 0;


920  for (i = 0; i < l; i++)


921  sum_alpha += alpha[i];


922 


923  if (Cp == Cn)


924  Procedures.info("nu = " + sum_alpha / (Cp * prob.Count) + "\n");


925 


926  for (i = 0; i < l; i++)


927  alpha[i] *= y[i];


928  }


929 


930  private static void solve_nu_svc(Problem prob, Parameter param,


931  double[] alpha, Solver.SolutionInfo si) {


932  int i;


933  int l = prob.Count;


934  double nu = param.Nu;


935 


936  sbyte[] y = new sbyte[l];


937 


938  for (i = 0; i < l; i++)


939  if (prob.Y[i] > 0)


940  y[i] = +1;


941  else


942  y[i] = 1;


943 


944  double sum_pos = nu * l / 2;


945  double sum_neg = nu * l / 2;


946 


947  for (i = 0; i < l; i++)


948  if (y[i] == +1) {


949  alpha[i] = Math.Min(1.0, sum_pos);


950  sum_pos = alpha[i];


951  } else {


952  alpha[i] = Math.Min(1.0, sum_neg);


953  sum_neg = alpha[i];


954  }


955 


956  double[] zeros = new double[l];


957 


958  for (i = 0; i < l; i++)


959  zeros[i] = 0;


960 


961  Solver_NU s = new Solver_NU();


962  s.Solve(l, new SVC_Q(prob, param, y), zeros, y, alpha, 1.0, 1.0, param.EPS, si, param.Shrinking);


963  double r = si.r;


964 


965  Procedures.info("C = " + 1 / r + "\n");


966 


967  for (i = 0; i < l; i++)


968  alpha[i] *= y[i] / r;


969 


970  si.rho /= r;


971  si.obj /= (r * r);


972  si.upper_bound_p = 1 / r;


973  si.upper_bound_n = 1 / r;


974  }


975 


976  private static void solve_one_class(Problem prob, Parameter param,


977  double[] alpha, Solver.SolutionInfo si) {


978  int l = prob.Count;


979  double[] zeros = new double[l];


980  sbyte[] ones = new sbyte[l];


981  int i;


982 


983  int n = (int)(param.Nu * prob.Count); // # of alpha's at upper bound


984 


985  for (i = 0; i < n; i++)


986  alpha[i] = 1;


987  if (n < prob.Count)


988  alpha[n] = param.Nu * prob.Count  n;


989  for (i = n + 1; i < l; i++)


990  alpha[i] = 0;


991 


992  for (i = 0; i < l; i++) {


993  zeros[i] = 0;


994  ones[i] = 1;


995  }


996 


997  Solver s = new Solver();


998  s.Solve(l, new ONE_CLASS_Q(prob, param), zeros, ones, alpha, 1.0, 1.0, param.EPS, si, param.Shrinking);


999  }


1000 


1001  private static void solve_epsilon_svr(Problem prob, Parameter param, double[] alpha, Solver.SolutionInfo si) {


1002  int l = prob.Count;


1003  double[] alpha2 = new double[2 * l];


1004  double[] linear_term = new double[2 * l];


1005  sbyte[] y = new sbyte[2 * l];


1006  int i;


1007 


1008  for (i = 0; i < l; i++) {


1009  alpha2[i] = 0;


1010  linear_term[i] = param.P  prob.Y[i];


1011  y[i] = 1;


1012 


1013  alpha2[i + l] = 0;


1014  linear_term[i + l] = param.P + prob.Y[i];


1015  y[i + l] = 1;


1016  }


1017 


1018  Solver s = new Solver();


1019  s.Solve(2 * l, new SVR_Q(prob, param), linear_term, y, alpha2, param.C, param.C, param.EPS, si, param.Shrinking);


1020 


1021  double sum_alpha = 0;


1022  for (i = 0; i < l; i++) {


1023  alpha[i] = alpha2[i]  alpha2[i + l];


1024  sum_alpha += Math.Abs(alpha[i]);


1025  }


1026  Procedures.info("nu = " + sum_alpha / (param.C * l) + "\n");


1027  }


1028 


1029  private static void solve_nu_svr(Problem prob, Parameter param,


1030  double[] alpha, Solver.SolutionInfo si) {


1031  int l = prob.Count;


1032  double C = param.C;


1033  double[] alpha2 = new double[2 * l];


1034  double[] linear_term = new double[2 * l];


1035  sbyte[] y = new sbyte[2 * l];


1036  int i;


1037 


1038  double sum = C * param.Nu * l / 2;


1039  for (i = 0; i < l; i++) {


1040  alpha2[i] = alpha2[i + l] = Math.Min(sum, C);


1041  sum = alpha2[i];


1042 


1043  linear_term[i] = prob.Y[i];


1044  y[i] = 1;


1045 


1046  linear_term[i + l] = prob.Y[i];


1047  y[i + l] = 1;


1048  }


1049 


1050  Solver_NU s = new Solver_NU();


1051  s.Solve(2 * l, new SVR_Q(prob, param), linear_term, y, alpha2, C, C, param.EPS, si, param.Shrinking);


1052 


1053  Procedures.info("epsilon = " + (si.r) + "\n");


1054 


1055  for (i = 0; i < l; i++)


1056  alpha[i] = alpha2[i]  alpha2[i + l];


1057  }


1058 


1059  //


1060  // decision_function


1061  //


1062  internal class decision_function {


1063  public double[] alpha;


1064  public double rho;


1065  };


1066 


1067  static decision_function svm_train_one(Problem prob, Parameter param, double Cp, double Cn) {


1068  double[] alpha = new double[prob.Count];


1069  Solver.SolutionInfo si = new Solver.SolutionInfo();


1070  switch (param.SvmType) {


1071  case SvmType.C_SVC:


1072  solve_c_svc(prob, param, alpha, si, Cp, Cn);


1073  break;


1074  case SvmType.NU_SVC:


1075  solve_nu_svc(prob, param, alpha, si);


1076  break;


1077  case SvmType.ONE_CLASS:


1078  solve_one_class(prob, param, alpha, si);


1079  break;


1080  case SvmType.EPSILON_SVR:


1081  solve_epsilon_svr(prob, param, alpha, si);


1082  break;


1083  case SvmType.NU_SVR:


1084  solve_nu_svr(prob, param, alpha, si);


1085  break;


1086  }


1087 


1088  Procedures.info("obj = " + si.obj + ", rho = " + si.rho + "\n");


1089 


1090  // output SVs


1091 


1092  int nSV = 0;


1093  int nBSV = 0;


1094  for (int i = 0; i < prob.Count; i++) {


1095  if (Math.Abs(alpha[i]) > 0) {


1096  ++nSV;


1097  if (prob.Y[i] > 0) {


1098  if (Math.Abs(alpha[i]) >= si.upper_bound_p)


1099  ++nBSV;


1100  } else {


1101  if (Math.Abs(alpha[i]) >= si.upper_bound_n)


1102  ++nBSV;


1103  }


1104  }


1105  }


1106 


1107  Procedures.info("nSV = " + nSV + ", nBSV = " + nBSV + "\n");


1108 


1109  decision_function f = new decision_function();


1110  f.alpha = alpha;


1111  f.rho = si.rho;


1112  return f;


1113  }


1114 


1115  // Platt's binary SVM Probablistic Output: an improvement from Lin et al.


1116  private static void sigmoid_train(int l, double[] dec_values, double[] labels,


1117  double[] probAB) {


1118  double A, B;


1119  double prior1 = 0, prior0 = 0;


1120  int i;


1121 


1122  for (i = 0; i < l; i++)


1123  if (labels[i] > 0) prior1 += 1;


1124  else prior0 += 1;


1125 


1126  int Max_iter = 100; // Maximal number of iterations


1127  double Min_step = 1e10; // Minimal step taken in line search


1128  double sigma = 1e12; // For numerically strict PD of Hessian


1129  double eps = 1e5;


1130  double hiTarget = (prior1 + 1.0) / (prior1 + 2.0);


1131  double loTarget = 1 / (prior0 + 2.0);


1132  double[] t = new double[l];


1133  double fApB, p, q, h11, h22, h21, g1, g2, det, dA, dB, gd, stepsize;


1134  double newA, newB, newf, d1, d2;


1135  int iter;


1136 


1137  // Initial Point and Initial Fun Value


1138  A = 0.0; B = Math.Log((prior0 + 1.0) / (prior1 + 1.0));


1139  double fval = 0.0;


1140 


1141  for (i = 0; i < l; i++) {


1142  if (labels[i] > 0) t[i] = hiTarget;


1143  else t[i] = loTarget;


1144  fApB = dec_values[i] * A + B;


1145  if (fApB >= 0)


1146  fval += t[i] * fApB + Math.Log(1 + Math.Exp(fApB));


1147  else


1148  fval += (t[i]  1) * fApB + Math.Log(1 + Math.Exp(fApB));


1149  }


1150  for (iter = 0; iter < Max_iter; iter++) {


1151  // Update Gradient and Hessian (use H' = H + sigma I)


1152  h11 = sigma; // numerically ensures strict PD


1153  h22 = sigma;


1154  h21 = 0.0; g1 = 0.0; g2 = 0.0;


1155  for (i = 0; i < l; i++) {


1156  fApB = dec_values[i] * A + B;


1157  if (fApB >= 0) {


1158  p = Math.Exp(fApB) / (1.0 + Math.Exp(fApB));


1159  q = 1.0 / (1.0 + Math.Exp(fApB));


1160  } else {


1161  p = 1.0 / (1.0 + Math.Exp(fApB));


1162  q = Math.Exp(fApB) / (1.0 + Math.Exp(fApB));


1163  }


1164  d2 = p * q;


1165  h11 += dec_values[i] * dec_values[i] * d2;


1166  h22 += d2;


1167  h21 += dec_values[i] * d2;


1168  d1 = t[i]  p;


1169  g1 += dec_values[i] * d1;


1170  g2 += d1;


1171  }


1172 


1173  // Stopping Criteria


1174  if (Math.Abs(g1) < eps && Math.Abs(g2) < eps)


1175  break;


1176 


1177  // Finding Newton direction: inv(H') * g


1178  det = h11 * h22  h21 * h21;


1179  dA = (h22 * g1  h21 * g2) / det;


1180  dB = (h21 * g1 + h11 * g2) / det;


1181  gd = g1 * dA + g2 * dB;


1182 


1183 


1184  stepsize = 1; // Line Search


1185  while (stepsize >= Min_step) {


1186  newA = A + stepsize * dA;


1187  newB = B + stepsize * dB;


1188 


1189  // New function value


1190  newf = 0.0;


1191  for (i = 0; i < l; i++) {


1192  fApB = dec_values[i] * newA + newB;


1193  if (fApB >= 0)


1194  newf += t[i] * fApB + Math.Log(1 + Math.Exp(fApB));


1195  else


1196  newf += (t[i]  1) * fApB + Math.Log(1 + Math.Exp(fApB));


1197  }


1198  // Check sufficient decrease


1199  if (newf < fval + 0.0001 * stepsize * gd) {


1200  A = newA; B = newB; fval = newf;


1201  break;


1202  } else


1203  stepsize = stepsize / 2.0;


1204  }


1205 


1206  if (stepsize < Min_step) {


1207  Procedures.info("Line search fails in twoclass probability estimates\n");


1208  break;


1209  }


1210  }


1211 


1212  if (iter >= Max_iter)


1213  Procedures.info("Reaching Maximal iterations in twoclass probability estimates\n");


1214  probAB[0] = A; probAB[1] = B;


1215  }


1216 


1217  private static double sigmoid_predict(double decision_value, double A, double B) {


1218  double fApB = decision_value * A + B;


1219  if (fApB >= 0)


1220  return Math.Exp(fApB) / (1.0 + Math.Exp(fApB));


1221  else


1222  return 1.0 / (1 + Math.Exp(fApB));


1223  }


1224 


1225  // Method 2 from the multiclass_prob paper by Wu, Lin, and Weng


1226  private static void multiclass_probability(int k, double[,] r, double[] p) {


1227  int t, j;


1228  int iter = 0, Max_iter = Math.Max(100, k);


1229  double[,] Q = new double[k, k];


1230  double[] Qp = new double[k];


1231  double pQp, eps = 0.005 / k;


1232 


1233  for (t = 0; t < k; t++) {


1234  p[t] = 1.0 / k; // Valid if k = 1


1235  Q[t, t] = 0;


1236  for (j = 0; j < t; j++) {


1237  Q[t, t] += r[j, t] * r[j, t];


1238  Q[t, j] = Q[j, t];


1239  }


1240  for (j = t + 1; j < k; j++) {


1241  Q[t, t] += r[j, t] * r[j, t];


1242  Q[t, j] = r[j, t] * r[t, j];


1243  }


1244  }


1245  for (iter = 0; iter < Max_iter; iter++) {


1246  // stopping condition, recalculate QP,pQP for numerical accuracy


1247  pQp = 0;


1248  for (t = 0; t < k; t++) {


1249  Qp[t] = 0;


1250  for (j = 0; j < k; j++)


1251  Qp[t] += Q[t, j] * p[j];


1252  pQp += p[t] * Qp[t];


1253  }


1254  double Max_error = 0;


1255  for (t = 0; t < k; t++) {


1256  double error = Math.Abs(Qp[t]  pQp);


1257  if (error > Max_error)


1258  Max_error = error;


1259  }


1260  if (Max_error < eps) break;


1261 


1262  for (t = 0; t < k; t++) {


1263  double diff = (Qp[t] + pQp) / Q[t, t];


1264  p[t] += diff;


1265  pQp = (pQp + diff * (diff * Q[t, t] + 2 * Qp[t])) / (1 + diff) / (1 + diff);


1266  for (j = 0; j < k; j++) {


1267  Qp[j] = (Qp[j] + diff * Q[t, j]) / (1 + diff);


1268  p[j] /= (1 + diff);


1269  }


1270  }


1271  }


1272  if (iter >= Max_iter)


1273  Procedures.info("Exceeds Max_iter in multiclass_prob\n");


1274  }


1275 


1276  // Crossvalidation decision values for probability estimates


1277  private static void svm_binary_svc_probability(Problem prob, Parameter param, double Cp, double Cn, double[] probAB) {


1278  int i;


1279  int nr_fold = 5;


1280  int[] perm = new int[prob.Count];


1281  double[] dec_values = new double[prob.Count];


1282 


1283  // random shuffle


1284  Random rand = new Random();


1285  for (i = 0; i < prob.Count; i++) perm[i] = i;


1286  for (i = 0; i < prob.Count; i++) {


1287  int j = i + (int)(rand.NextDouble() * (prob.Count  i));


1288  do { int _ = perm[i]; perm[i] = perm[j]; perm[j] = _; } while (false);


1289  }


1290  for (i = 0; i < nr_fold; i++) {


1291  int begin = i * prob.Count / nr_fold;


1292  int end = (i + 1) * prob.Count / nr_fold;


1293  int j, k;


1294  Problem subprob = new Problem();


1295 


1296  subprob.Count = prob.Count  (end  begin);


1297  subprob.X = new Node[subprob.Count][];


1298  subprob.Y = new double[subprob.Count];


1299 


1300  k = 0;


1301  for (j = 0; j < begin; j++) {


1302  subprob.X[k] = prob.X[perm[j]];


1303  subprob.Y[k] = prob.Y[perm[j]];


1304  ++k;


1305  }


1306  for (j = end; j < prob.Count; j++) {


1307  subprob.X[k] = prob.X[perm[j]];


1308  subprob.Y[k] = prob.Y[perm[j]];


1309  ++k;


1310  }


1311  int p_count = 0, n_count = 0;


1312  for (j = 0; j < k; j++)


1313  if (subprob.Y[j] > 0)


1314  p_count++;


1315  else


1316  n_count++;


1317 


1318  if (p_count == 0 && n_count == 0)


1319  for (j = begin; j < end; j++)


1320  dec_values[perm[j]] = 0;


1321  else if (p_count > 0 && n_count == 0)


1322  for (j = begin; j < end; j++)


1323  dec_values[perm[j]] = 1;


1324  else if (p_count == 0 && n_count > 0)


1325  for (j = begin; j < end; j++)


1326  dec_values[perm[j]] = 1;


1327  else {


1328  Parameter subparam = (Parameter)param.Clone();


1329  subparam.Probability = false;


1330  subparam.C = 1.0;


1331  subparam.Weights[1] = Cp;


1332  subparam.Weights[1] = Cn;


1333  Model submodel = svm_train(subprob, subparam);


1334  for (j = begin; j < end; j++) {


1335  double[] dec_value = new double[1];


1336  svm_predict_values(submodel, prob.X[perm[j]], dec_value);


1337  dec_values[perm[j]] = dec_value[0];


1338  // ensure +1 1 order; reason not using CV subroutine


1339  dec_values[perm[j]] *= submodel.ClassLabels[0];


1340  }


1341  }


1342  }


1343  sigmoid_train(prob.Count, dec_values, prob.Y, probAB);


1344  }


1345 


1346  // Return parameter of a Laplace distribution


1347  private static double svm_svr_probability(Problem prob, Parameter param) {


1348  int i;


1349  int nr_fold = 5;


1350  double[] ymv = new double[prob.Count];


1351  double mae = 0;


1352 


1353  Parameter newparam = (Parameter)param.Clone();


1354  newparam.Probability = false;


1355  svm_cross_validation(prob, newparam, nr_fold, ymv, true);


1356  for (i = 0; i < prob.Count; i++) {


1357  ymv[i] = prob.Y[i]  ymv[i];


1358  mae += Math.Abs(ymv[i]);


1359  }


1360  mae /= prob.Count;


1361  double std = Math.Sqrt(2 * mae * mae);


1362  int count = 0;


1363  mae = 0;


1364  for (i = 0; i < prob.Count; i++)


1365  if (Math.Abs(ymv[i]) > 5 * std)


1366  count = count + 1;


1367  else


1368  mae += Math.Abs(ymv[i]);


1369  mae /= (prob.Count  count);


1370  Procedures.info("Prob. model for test data: target value = predicted value + z,\nz: Laplace distribution e^(z/sigma)/(2sigma),sigma=" + mae + "\n");


1371  return mae;


1372  }


1373 


1374  // label: label name, start: begin of each class, count: #data of classes, perm: indices to the original data


1375  // perm, length l, must be allocated before calling this subroutine


1376  private static void svm_group_classes(Problem prob, int[] nr_class_ret, int[][] label_ret, int[][] start_ret, int[][] count_ret, int[] perm) {


1377  int l = prob.Count;


1378  int Max_nr_class = 16;


1379  int nr_class = 0;


1380  int[] label = new int[Max_nr_class];


1381  int[] count = new int[Max_nr_class];


1382  int[] data_label = new int[l];


1383  int i;


1384 


1385  for (i = 0; i < l; i++) {


1386  int this_label = (int)(prob.Y[i]);


1387  int j;


1388  for (j = 0; j < nr_class; j++) {


1389  if (this_label == label[j]) {


1390  ++count[j];


1391  break;


1392  }


1393  }


1394  data_label[i] = j;


1395  if (j == nr_class) {


1396  if (nr_class == Max_nr_class) {


1397  Max_nr_class *= 2;


1398  int[] new_data = new int[Max_nr_class];


1399  Array.Copy(label, 0, new_data, 0, label.Length);


1400  label = new_data;


1401  new_data = new int[Max_nr_class];


1402  Array.Copy(count, 0, new_data, 0, count.Length);


1403  count = new_data;


1404  }


1405  label[nr_class] = this_label;


1406  count[nr_class] = 1;


1407  ++nr_class;


1408  }


1409  }


1410 


1411  int[] start = new int[nr_class];


1412  start[0] = 0;


1413  for (i = 1; i < nr_class; i++)


1414  start[i] = start[i  1] + count[i  1];


1415  for (i = 0; i < l; i++) {


1416  perm[start[data_label[i]]] = i;


1417  ++start[data_label[i]];


1418  }


1419  start[0] = 0;


1420  for (i = 1; i < nr_class; i++)


1421  start[i] = start[i  1] + count[i  1];


1422 


1423  nr_class_ret[0] = nr_class;


1424  label_ret[0] = label;


1425  start_ret[0] = start;


1426  count_ret[0] = count;


1427  }


1428 


1429  //


1430  // Interface functions


1431  //


1432  public static Model svm_train(Problem prob, Parameter param) {


1433  Model model = new Model();


1434  model.Parameter = param;


1435 


1436  if (param.SvmType == SvmType.ONE_CLASS 


1437  param.SvmType == SvmType.EPSILON_SVR 


1438  param.SvmType == SvmType.NU_SVR) {


1439  // regression or oneclasssvm


1440  model.NumberOfClasses = 2;


1441  model.ClassLabels = null;


1442  model.NumberOfSVPerClass = null;


1443  model.PairwiseProbabilityA = null; model.PairwiseProbabilityB = null;


1444  model.SupportVectorCoefficients = new double[1][];


1445 


1446  if (param.Probability &&


1447  (param.SvmType == SvmType.EPSILON_SVR 


1448  param.SvmType == SvmType.NU_SVR)) {


1449  model.PairwiseProbabilityA = new double[1];


1450  model.PairwiseProbabilityA[0] = svm_svr_probability(prob, param);


1451  }


1452 


1453  decision_function f = svm_train_one(prob, param, 0, 0);


1454  model.Rho = new double[1];


1455  model.Rho[0] = f.rho;


1456 


1457  int nSV = 0;


1458  int i;


1459  for (i = 0; i < prob.Count; i++)


1460  if (Math.Abs(f.alpha[i]) > 0) ++nSV;


1461  model.SupportVectorCount = nSV;


1462  model.SupportVectorIndizes = new int[nSV];


1463  model.SupportVectors = new Node[nSV][];


1464  model.SupportVectorCoefficients[0] = new double[nSV];


1465 


1466  int j = 0;


1467  for (i = 0; i < prob.Count; i++)


1468  if (Math.Abs(f.alpha[i]) > 0) {


1469  model.SupportVectors[j] = prob.X[i];


1470  model.SupportVectorIndizes[j] = i;


1471  model.SupportVectorCoefficients[0][j] = f.alpha[i];


1472 


1473  ++j;


1474  }


1475  } else {


1476  // classification


1477  int l = prob.Count;


1478  int[] tmp_nr_class = new int[1];


1479  int[][] tmp_label = new int[1][];


1480  int[][] tmp_start = new int[1][];


1481  int[][] tmp_count = new int[1][];


1482  int[] perm = new int[l];


1483 


1484  // group training data of the same class


1485  svm_group_classes(prob, tmp_nr_class, tmp_label, tmp_start, tmp_count, perm);


1486  int nr_class = tmp_nr_class[0];


1487  int[] label = tmp_label[0];


1488  int[] start = tmp_start[0];


1489  int[] count = tmp_count[0];


1490  Node[][] x = new Node[l][];


1491  int i;


1492  for (i = 0; i < l; i++)


1493  x[i] = prob.X[perm[i]];


1494 


1495  // calculate weighted C


1496 


1497  double[] weighted_C = new double[nr_class];


1498  for (i = 0; i < nr_class; i++)


1499  weighted_C[i] = param.C;


1500  foreach (int weightedLabel in param.Weights.Keys) {


1501  int index = Array.IndexOf<int>(label, weightedLabel);


1502  if (index < 0)


1503  Console.Error.WriteLine("warning: class label " + weightedLabel + " specified in weight is not found");


1504  else weighted_C[index] *= param.Weights[weightedLabel];


1505  }


1506 


1507  // train k*(k1)/2 models


1508 


1509  bool[] nonzero = new bool[l];


1510  for (i = 0; i < l; i++)


1511  nonzero[i] = false;


1512  decision_function[] f = new decision_function[nr_class * (nr_class  1) / 2];


1513 


1514  double[] probA = null, probB = null;


1515  if (param.Probability) {


1516  probA = new double[nr_class * (nr_class  1) / 2];


1517  probB = new double[nr_class * (nr_class  1) / 2];


1518  }


1519 


1520  int p = 0;


1521  for (i = 0; i < nr_class; i++)


1522  for (int j = i + 1; j < nr_class; j++) {


1523  Problem sub_prob = new Problem();


1524  int si = start[i], sj = start[j];


1525  int ci = count[i], cj = count[j];


1526  sub_prob.Count = ci + cj;


1527  sub_prob.X = new Node[sub_prob.Count][];


1528  sub_prob.Y = new double[sub_prob.Count];


1529  int k;


1530  for (k = 0; k < ci; k++) {


1531  sub_prob.X[k] = x[si + k];


1532  sub_prob.Y[k] = +1;


1533  }


1534  for (k = 0; k < cj; k++) {


1535  sub_prob.X[ci + k] = x[sj + k];


1536  sub_prob.Y[ci + k] = 1;


1537  }


1538 


1539  if (param.Probability) {


1540  double[] probAB = new double[2];


1541  svm_binary_svc_probability(sub_prob, param, weighted_C[i], weighted_C[j], probAB);


1542  probA[p] = probAB[0];


1543  probB[p] = probAB[1];


1544  }


1545 


1546  f[p] = svm_train_one(sub_prob, param, weighted_C[i], weighted_C[j]);


1547  for (k = 0; k < ci; k++)


1548  if (!nonzero[si + k] && Math.Abs(f[p].alpha[k]) > 0)


1549  nonzero[si + k] = true;


1550  for (k = 0; k < cj; k++)


1551  if (!nonzero[sj + k] && Math.Abs(f[p].alpha[ci + k]) > 0)


1552  nonzero[sj + k] = true;


1553  ++p;


1554  }


1555 


1556  // build output


1557 


1558  model.NumberOfClasses = nr_class;


1559 


1560  model.ClassLabels = new int[nr_class];


1561  for (i = 0; i < nr_class; i++)


1562  model.ClassLabels[i] = label[i];


1563 


1564  model.Rho = new double[nr_class * (nr_class  1) / 2];


1565  for (i = 0; i < nr_class * (nr_class  1) / 2; i++)


1566  model.Rho[i] = f[i].rho;


1567 


1568  if (param.Probability) {


1569  model.PairwiseProbabilityA = new double[nr_class * (nr_class  1) / 2];


1570  model.PairwiseProbabilityB = new double[nr_class * (nr_class  1) / 2];


1571  for (i = 0; i < nr_class * (nr_class  1) / 2; i++) {


1572  model.PairwiseProbabilityA[i] = probA[i];


1573  model.PairwiseProbabilityB[i] = probB[i];


1574  }


1575  } else {


1576  model.PairwiseProbabilityA = null;


1577  model.PairwiseProbabilityB = null;


1578  }


1579 


1580  int nnz = 0;


1581  int[] nz_count = new int[nr_class];


1582  model.NumberOfSVPerClass = new int[nr_class];


1583  for (i = 0; i < nr_class; i++) {


1584  int nSV = 0;


1585  for (int j = 0; j < count[i]; j++)


1586  if (nonzero[start[i] + j]) {


1587  ++nSV;


1588  ++nnz;


1589  }


1590  model.NumberOfSVPerClass[i] = nSV;


1591  nz_count[i] = nSV;


1592  }


1593 


1594  Procedures.info("Total nSV = " + nnz + "\n");


1595 


1596  model.SupportVectorCount = nnz;


1597  model.SupportVectors = new Node[nnz][];


1598  model.SupportVectorIndizes = new int[nnz];


1599  p = 0;


1600  for (i = 0; i < l; i++) {


1601  if (nonzero[i]) {


1602  model.SupportVectors[p] = x[i];


1603  model.SupportVectorIndizes[p] = i;


1604  p++;


1605  }


1606  }


1607 


1608  int[] nz_start = new int[nr_class];


1609  nz_start[0] = 0;


1610  for (i = 1; i < nr_class; i++)


1611  nz_start[i] = nz_start[i  1] + nz_count[i  1];


1612 


1613  model.SupportVectorCoefficients = new double[nr_class  1][];


1614  for (i = 0; i < nr_class  1; i++)


1615  model.SupportVectorCoefficients[i] = new double[nnz];


1616 


1617  p = 0;


1618  for (i = 0; i < nr_class; i++)


1619  for (int j = i + 1; j < nr_class; j++) {


1620  // classifier (i,j): coefficients with


1621  // i are in sv_coef[j1][nz_start[i]...],


1622  // j are in sv_coef[i][nz_start[j]...]


1623 


1624  int si = start[i];


1625  int sj = start[j];


1626  int ci = count[i];


1627  int cj = count[j];


1628 


1629  int q = nz_start[i];


1630  int k;


1631  for (k = 0; k < ci; k++)


1632  if (nonzero[si + k])


1633  model.SupportVectorCoefficients[j  1][q++] = f[p].alpha[k];


1634  q = nz_start[j];


1635  for (k = 0; k < cj; k++)


1636  if (nonzero[sj + k])


1637  model.SupportVectorCoefficients[i][q++] = f[p].alpha[ci + k];


1638  ++p;


1639  }


1640  }


1641  return model;


1642  }


1643 


1644  // Stratified cross validation


1645  public static void svm_cross_validation(Problem prob, Parameter param, int nr_fold, double[] target, bool shuffleTraining) {


1646  Random rand = new Random();


1647  int i;


1648  int[] fold_start = new int[nr_fold + 1];


1649  int l = prob.Count;


1650  int[] perm = new int[l];


1651 


1652  // stratified cv may not give leaveoneout rate


1653  // Each class to l folds > some folds may have zero elements


1654  if ((param.SvmType == SvmType.C_SVC 


1655  param.SvmType == SvmType.NU_SVC) && nr_fold < l) {


1656  int[] tmp_nr_class = new int[1];


1657  int[][] tmp_label = new int[1][];


1658  int[][] tmp_start = new int[1][];


1659  int[][] tmp_count = new int[1][];


1660 


1661  svm_group_classes(prob, tmp_nr_class, tmp_label, tmp_start, tmp_count, perm);


1662 


1663  int nr_class = tmp_nr_class[0];


1664  int[] label = tmp_label[0];


1665  int[] start = tmp_start[0];


1666  int[] count = tmp_count[0];


1667 


1668  // random shuffle and then data grouped by fold using the array perm


1669  int[] fold_count = new int[nr_fold];


1670  int c;


1671  int[] index = new int[l];


1672  for (i = 0; i < l; i++)


1673  index[i] = perm[i];


1674  for (c = 0; c < nr_class; c++)


1675  for (i = 0; i < count[c]; i++) {


1676  int j = i + (int)(rand.NextDouble() * (count[c]  i));


1677  do { int _ = index[start[c] + j]; index[start[c] + j] = index[start[c] + i]; index[start[c] + i] = _; } while (false);


1678  }


1679  for (i = 0; i < nr_fold; i++) {


1680  fold_count[i] = 0;


1681  for (c = 0; c < nr_class; c++)


1682  fold_count[i] += (i + 1) * count[c] / nr_fold  i * count[c] / nr_fold;


1683  }


1684  fold_start[0] = 0;


1685  for (i = 1; i <= nr_fold; i++)


1686  fold_start[i] = fold_start[i  1] + fold_count[i  1];


1687  for (c = 0; c < nr_class; c++)


1688  for (i = 0; i < nr_fold; i++) {


1689  int begin = start[c] + i * count[c] / nr_fold;


1690  int end = start[c] + (i + 1) * count[c] / nr_fold;


1691  for (int j = begin; j < end; j++) {


1692  perm[fold_start[i]] = index[j];


1693  fold_start[i]++;


1694  }


1695  }


1696  fold_start[0] = 0;


1697  for (i = 1; i <= nr_fold; i++)


1698  fold_start[i] = fold_start[i  1] + fold_count[i  1];


1699  } else {


1700  for (i = 0; i < l; i++) perm[i] = i;


1701  if (shuffleTraining) {


1702  for (i = 0; i < l; i++) {


1703  int j = i + (int)(rand.NextDouble() * (l  i));


1704  do { int _ = perm[i]; perm[i] = perm[j]; perm[j] = _; } while (false);


1705  }


1706  }


1707  for (i = 0; i <= nr_fold; i++)


1708  fold_start[i] = i * l / nr_fold;


1709  }


1710 


1711  for (i = 0; i < nr_fold; i++) {


1712  int begin = fold_start[i];


1713  int end = fold_start[i + 1];


1714  int j, k;


1715  Problem subprob = new Problem();


1716 


1717  subprob.Count = l  (end  begin);


1718  subprob.X = new Node[subprob.Count][];


1719  subprob.Y = new double[subprob.Count];


1720 


1721  k = 0;


1722  for (j = 0; j < begin; j++) {


1723  subprob.X[k] = prob.X[perm[j]];


1724  subprob.Y[k] = prob.Y[perm[j]];


1725  ++k;


1726  }


1727  for (j = end; j < l; j++) {


1728  subprob.X[k] = prob.X[perm[j]];


1729  subprob.Y[k] = prob.Y[perm[j]];


1730  ++k;


1731  }


1732  Model submodel = svm_train(subprob, param);


1733  if (param.Probability &&


1734  (param.SvmType == SvmType.C_SVC 


1735  param.SvmType == SvmType.NU_SVC)) {


1736  double[] prob_estimates = new double[svm_get_nr_class(submodel)];


1737  for (j = begin; j < end; j++)


1738  target[perm[j]] = svm_predict_probability(submodel, prob.X[perm[j]], prob_estimates);


1739  } else


1740  for (j = begin; j < end; j++)


1741  target[perm[j]] = svm_predict(submodel, prob.X[perm[j]]);


1742  }


1743  }


1744 


1745  public static SvmType svm_get_svm_type(Model model) {


1746  return model.Parameter.SvmType;


1747  }


1748 


1749  public static int svm_get_nr_class(Model model) {


1750  return model.NumberOfClasses;


1751  }


1752 


1753  public static void svm_get_labels(Model model, int[] label) {


1754  if (model.ClassLabels != null)


1755  for (int i = 0; i < model.NumberOfClasses; i++)


1756  label[i] = model.ClassLabels[i];


1757  }


1758 


1759  public static double svm_get_svr_probability(Model model) {


1760  if ((model.Parameter.SvmType == SvmType.EPSILON_SVR  model.Parameter.SvmType == SvmType.NU_SVR) &&


1761  model.PairwiseProbabilityA != null)


1762  return model.PairwiseProbabilityA[0];


1763  else {


1764  Console.Error.WriteLine("Model doesn't contain information for SVR probability inference");


1765  return 0;


1766  }


1767  }


1768 


1769  public static void svm_predict_values(Model model, Node[] x, double[] dec_values) {


1770  if (model.Parameter.SvmType == SvmType.ONE_CLASS 


1771  model.Parameter.SvmType == SvmType.EPSILON_SVR 


1772  model.Parameter.SvmType == SvmType.NU_SVR) {


1773  double[] sv_coef = model.SupportVectorCoefficients[0];


1774  double sum = 0;


1775  for (int i = 0; i < model.SupportVectorCount; i++)


1776  sum += sv_coef[i] * Kernel.KernelFunction(x, model.SupportVectors[i], model.Parameter);


1777  sum = model.Rho[0];


1778  dec_values[0] = sum;


1779  } else {


1780  int i;


1781  int nr_class = model.NumberOfClasses;


1782  int l = model.SupportVectorCount;


1783 


1784  double[] kvalue = new double[l];


1785  for (i = 0; i < l; i++)


1786  kvalue[i] = Kernel.KernelFunction(x, model.SupportVectors[i], model.Parameter);


1787 


1788  int[] start = new int[nr_class];


1789  start[0] = 0;


1790  for (i = 1; i < nr_class; i++)


1791  start[i] = start[i  1] + model.NumberOfSVPerClass[i  1];


1792 


1793  int p = 0;


1794  for (i = 0; i < nr_class; i++)


1795  for (int j = i + 1; j < nr_class; j++) {


1796  double sum = 0;


1797  int si = start[i];


1798  int sj = start[j];


1799  int ci = model.NumberOfSVPerClass[i];


1800  int cj = model.NumberOfSVPerClass[j];


1801 


1802  int k;


1803  double[] coef1 = model.SupportVectorCoefficients[j  1];


1804  double[] coef2 = model.SupportVectorCoefficients[i];


1805  for (k = 0; k < ci; k++)


1806  sum += coef1[si + k] * kvalue[si + k];


1807  for (k = 0; k < cj; k++)


1808  sum += coef2[sj + k] * kvalue[sj + k];


1809  sum = model.Rho[p];


1810  dec_values[p] = sum;


1811  p++;


1812  }


1813  }


1814  }


1815 


1816  public static double svm_predict(Model model, Node[] x) {


1817  if (model.Parameter.SvmType == SvmType.ONE_CLASS 


1818  model.Parameter.SvmType == SvmType.EPSILON_SVR 


1819  model.Parameter.SvmType == SvmType.NU_SVR) {


1820  double[] res = new double[1];


1821  svm_predict_values(model, x, res);


1822 


1823  if (model.Parameter.SvmType == SvmType.ONE_CLASS)


1824  return (res[0] > 0) ? 1 : 1;


1825  else


1826  return res[0];


1827  } else {


1828  int i;


1829  int nr_class = model.NumberOfClasses;


1830  double[] dec_values = new double[nr_class * (nr_class  1) / 2];


1831  svm_predict_values(model, x, dec_values);


1832 


1833  int[] vote = new int[nr_class];


1834  for (i = 0; i < nr_class; i++)


1835  vote[i] = 0;


1836  int pos = 0;


1837  for (i = 0; i < nr_class; i++)


1838  for (int j = i + 1; j < nr_class; j++) {


1839  if (dec_values[pos++] > 0)


1840  ++vote[i];


1841  else


1842  ++vote[j];


1843  }


1844 


1845  int vote_Max_idx = 0;


1846  for (i = 1; i < nr_class; i++)


1847  if (vote[i] > vote[vote_Max_idx])


1848  vote_Max_idx = i;


1849  return model.ClassLabels[vote_Max_idx];


1850  }


1851  }


1852 


1853  public static double svm_predict_probability(Model model, Node[] x, double[] prob_estimates) {


1854  if ((model.Parameter.SvmType == SvmType.C_SVC  model.Parameter.SvmType == SvmType.NU_SVC) &&


1855  model.PairwiseProbabilityA != null && model.PairwiseProbabilityB != null) {


1856  int i;


1857  int nr_class = model.NumberOfClasses;


1858  double[] dec_values = new double[nr_class * (nr_class  1) / 2];


1859  svm_predict_values(model, x, dec_values);


1860 


1861  double Min_prob = 1e7;


1862  double[,] pairwise_prob = new double[nr_class, nr_class];


1863 


1864  int k = 0;


1865  for (i = 0; i < nr_class; i++) {


1866  for (int j = i + 1; j < nr_class; j++) {


1867  pairwise_prob[i, j] = Math.Min(Math.Max(sigmoid_predict(dec_values[k], model.PairwiseProbabilityA[k], model.PairwiseProbabilityB[k]), Min_prob), 1  Min_prob);


1868  pairwise_prob[j, i] = 1  pairwise_prob[i, j];


1869  k++;


1870  }


1871  }


1872  multiclass_probability(nr_class, pairwise_prob, prob_estimates);


1873 


1874  int prob_Max_idx = 0;


1875  for (i = 1; i < nr_class; i++)


1876  if (prob_estimates[i] > prob_estimates[prob_Max_idx])


1877  prob_Max_idx = i;


1878  return model.ClassLabels[prob_Max_idx];


1879  } else


1880  return svm_predict(model, x);


1881  }


1882 


1883  public static string svm_check_parameter(Problem prob, Parameter param) {


1884  // svm_type


1885 


1886  SvmType svm_type = param.SvmType;


1887 


1888  // kernel_type, degree


1889 


1890  KernelType kernel_type = param.KernelType;


1891 


1892  if (param.Degree < 0)


1893  return "degree of polynomial kernel < 0";


1894 


1895  // cache_size,eps,C,nu,p,shrinking


1896 


1897  if (param.CacheSize <= 0)


1898  return "cache_size <= 0";


1899 


1900  if (param.EPS <= 0)


1901  return "eps <= 0";


1902 


1903  if (param.Gamma == 0)


1904  param.Gamma = 1.0 / prob.MaxIndex;


1905 


1906  if (svm_type == SvmType.C_SVC 


1907  svm_type == SvmType.EPSILON_SVR 


1908  svm_type == SvmType.NU_SVR)


1909  if (param.C <= 0)


1910  return "C <= 0";


1911 


1912  if (svm_type == SvmType.NU_SVC 


1913  svm_type == SvmType.ONE_CLASS 


1914  svm_type == SvmType.NU_SVR)


1915  if (param.Nu <= 0  param.Nu > 1)


1916  return "nu <= 0 or nu > 1";


1917 


1918  if (svm_type == SvmType.EPSILON_SVR)


1919  if (param.P < 0)


1920  return "p < 0";


1921 


1922  if (param.Probability &&


1923  svm_type == SvmType.ONE_CLASS)


1924  return "oneclass SVM probability output not supported yet";


1925 


1926  // check whether nusvc is feasible


1927 


1928  if (svm_type == SvmType.NU_SVC) {


1929  int l = prob.Count;


1930  int Max_nr_class = 16;


1931  int nr_class = 0;


1932  int[] label = new int[Max_nr_class];


1933  int[] count = new int[Max_nr_class];


1934 


1935  int i;


1936  for (i = 0; i < l; i++) {


1937  int this_label = (int)prob.Y[i];


1938  int j;


1939  for (j = 0; j < nr_class; j++)


1940  if (this_label == label[j]) {


1941  ++count[j];


1942  break;


1943  }


1944 


1945  if (j == nr_class) {


1946  if (nr_class == Max_nr_class) {


1947  Max_nr_class *= 2;


1948  int[] new_data = new int[Max_nr_class];


1949  Array.Copy(label, 0, new_data, 0, label.Length);


1950  label = new_data;


1951 


1952  new_data = new int[Max_nr_class];


1953  Array.Copy(count, 0, new_data, 0, count.Length);


1954  count = new_data;


1955  }


1956  label[nr_class] = this_label;


1957  count[nr_class] = 1;


1958  ++nr_class;


1959  }


1960  }


1961 


1962  for (i = 0; i < nr_class; i++) {


1963  int n1 = count[i];


1964  for (int j = i + 1; j < nr_class; j++) {


1965  int n2 = count[j];


1966  if (param.Nu * (n1 + n2) / 2 > Math.Min(n1, n2))


1967  return "specified nu is infeasible";


1968  }


1969  }


1970  }


1971 


1972  return null;


1973  }


1974 


1975  public static int svm_check_probability_model(Model model) {


1976  if (((model.Parameter.SvmType == SvmType.C_SVC  model.Parameter.SvmType == SvmType.NU_SVC) &&


1977  model.PairwiseProbabilityA != null && model.PairwiseProbabilityB != null) 


1978  ((model.Parameter.SvmType == SvmType.EPSILON_SVR  model.Parameter.SvmType == SvmType.NU_SVR) &&


1979  model.PairwiseProbabilityA != null))


1980  return 1;


1981  else


1982  return 0;


1983  }


1984  }


1985  }

