smo.cc

Support Vector Machine smo

#line 825 "smo.w"

#line 947 "smo.w"

#include <vector>
#include <algorithm>
#include <functional>

#line 1475 "smo.w"
#include <cmath>

#line 1586 "smo.w"
#include <iostream>
#include <cstdlib>
#include <unistd.h>

#line 1710 "smo.w"
#include <string>
#include <vector>
#include <iostream>
#include <fstream>
#include <strstream>

#line 825 "smo.w"

using namespace std;

#line 901 "smo.w"

  int N = 0;                    /* N points(rows) */
  int d = -1;                   /* d variables */
  float C=0.05;
  float tolerance=0.001;
  float eps=0.001;
  float two_sigma_squared=2;

  vector<float> alph;           /* Lagrange multipliers */
  float b;                      /* threshold */
  vector<float> w;              /* weight vector: only for linear kernel */

  vector<float> error_cache;

  struct sparse_binary_vector {
      vector<int> id;
    };
  struct sparse_vector {
      vector<int> id;
      vector<float> val;
    };
  typedef vector<float> dense_vector;

  bool is_sparse_data = false;
  bool is_binary = false;
  /* use only one of these */
  vector<sparse_binary_vector> sparse_binary_points;
  vector<sparse_vector> sparse_points;
  vector<dense_vector> dense_points;

  vector<int> target;           /* class labels of training data points */
  bool is_test_only = false;
  bool is_linear_kernel = false;

  /* data points with index in [first_test_i .. N)
   * will be tested to compute error rate
   */
  int first_test_i = 0;

  /*
   * support vectors are within [0..end_support_i)
   */
  int end_support_i = -1;

#line 1015 "smo.w"
int takeStep(int i1, int i2);

#line 1135 "smo.w"
 float (*learned_func)(int) = NULL;

#line 1175 "smo.w"

 float (*kernel_func)(int,int)=NULL;

#line 1213 "smo.w"
 float delta_b;

#line 1351 "smo.w"

 float (*dot_product_func)(int,int)=NULL;

#line 1488 "smo.w"

    vector<float> precomputed_self_dot_product;

#line 827 "smo.w"


#line 963 "smo.w"
int examineExample(int i1)
{
  float y1, alph1, E1, r1;

  y1 = target[i1];
  alph1 = alph[i1];

  if (alph1 > 0 && alph1 < C)
     E1 = error_cache[i1];
  else 
     E1 = learned_func(i1) - y1;

  r1 = y1 * E1;
  if ((r1 < -tolerance && alph1 < C)
     || (r1 > tolerance && alph1 > 0))
    { 
      /* Try i2 by three ways; if successful, then immediately return 1; */
      
#line 991 "smo.w"
      {
        int k, i2;
        float tmax;
      
        for (i2 = (-1), tmax = 0, k = 0; k < end_support_i; k++)
          if (alph[k] > 0 && alph[k] < C) {
              float E2, temp;
      
              E2 = error_cache[k];
              temp = fabs(E1 - E2);
              if (temp > tmax)
                {
                  tmax = temp;
                  i2 = k;
                }
            }
      
        if (i2 >= 0) {
            if (takeStep (i1, i2))
              return 1;
          }
      }
#line 980 "smo.w"

      
#line 1021 "smo.w"
      {
        int k, k0;
        int i2;
      
        for (k0 = (int) (drand48 () * end_support_i), k = k0; k < end_support_i + k0; k++) {
            i2 = k % end_support_i;
            if (alph[i2] > 0 && alph[i2] < C) {
                if (takeStep(i1, i2))
                    return 1;
              }
          }
      }
#line 981 "smo.w"

      
#line 1038 "smo.w"
      {
        int k0, k, i2;
      
        for (k0 = (int)(drand48 () * end_support_i), k = k0; k < end_support_i + k0; k++) {
            i2 = k % end_support_i;
            if (takeStep(i1, i2))
                return 1;
          }
      }
#line 982 "smo.w"

    }

  return 0;
}
#line 1054 "smo.w"
int takeStep(int i1, int i2) { 
  int y1, y2, s;
  float alph1, alph2; /* old_values of alpha_1, alpha_2 */
  float a1, a2;       /* new values of alpha_1, alpha_2 */
  float E1, E2, L, H, k11, k22, k12, eta, Lobj, Hobj;

  if (i1 == i2) return 0;

  
#line 1120 "smo.w"
    alph1 = alph[i1];
    y1 = target[i1];
    if (alph1 > 0 && alph1 < C)
       E1 = error_cache[i1];
    else 
       E1 = learned_func(i1) - y1;
  
    alph2 = alph[i2];
    y2 = target[i2];
    if (alph2 > 0 && alph2 < C)
       E2 = error_cache[i2];
    else 
       E2 = learned_func(i2) - y2;
  
#line 1062 "smo.w"

  s = y1 * y2;

  
#line 1142 "smo.w"
  
     if (y1 == y2) {
       float gamma = alph1 + alph2;
       if (gamma > C) {
           L = gamma-C;
           H = C;
       }
       else {
           L = 0;
           H = gamma;
       }
     }
     else {
         float gamma = alph1 - alph2;
         if (gamma > 0) {
             L = 0;
             H = C - gamma;
         }
         else {
             L = -gamma;
             H = C;
         }
     }
  
#line 1065 "smo.w"

  if (L == H)
    return 0;

  
#line 1168 "smo.w"
    k11 = kernel_func(i1, i1);
    k12 = kernel_func(i1, i2);
    k22 = kernel_func(i2, i2);
    eta = 2 * k12 - k11 - k22;
  
#line 1069 "smo.w"

  if (eta < 0) {
      a2 = alph2 + y2 * (E2 - E1) / eta;
      if (a2 < L)
        a2 = L;
      else if (a2 > H)
        a2 = H;
    }
  else {
    
#line 1181 "smo.w"
    {
        float c1 = eta/2;
        float c2 = y2 * (E1-E2)- eta * alph2;
        Lobj = c1 * L * L + c2 * L;
        Hobj = c1 * H * H + c2 * H;
    }
#line 1078 "smo.w"

    if (Lobj > Hobj+eps)
      a2 = L;
    else if (Lobj < Hobj-eps)
      a2 = H;
    else
      a2 = alph2;
  }

  if (fabs(a2-alph2) < eps*(a2+alph2+eps))
    return 0;

  a1 = alph1 - s * (a2 - alph2);
  if (a1 < 0) {
        a2 += s * a1;
        a1 = 0;
  }
  else if (a1 > C) {
        float t = a1-C;
        a2 += s * t;
        a1 = C;
  }

  
#line 1193 "smo.w"
  {
    float b1, b2, bnew;
  
    if (a1 > 0 && a1 < C)
      bnew = b + E1 + y1 * (a1 - alph1) * k11 + y2 * (a2 - alph2) * k12;
    else {
        if (a2 > 0 && a2 < C)
            bnew = b + E2 + y1 * (a1 - alph1) * k12 + y2 * (a2 - alph2) * k22;
        else {
            b1 = b + E1 + y1 * (a1 - alph1) * k11 + y2 * (a2 - alph2) * k12;
            b2 = b + E2 + y1 * (a1 - alph1) * k12 + y2 * (a2 - alph2) * k22;
            bnew = (b1 + b2) / 2;
        }
    }
  
    delta_b = bnew - b;
    b = bnew;
  }
#line 1101 "smo.w"

  
#line 1226 "smo.w"
  
    if (is_linear_kernel) {
      float t1 = y1 * (a1 - alph1);
      float t2 = y2 * (a2 - alph2);
  
      if (is_sparse_data && is_binary) {
          int p1,num1,p2,num2;
  
          num1 = sparse_binary_points[i1].id.size();
          for (p1=0; p1<num1; p1++)
              w[sparse_binary_points[i1].id[p1]] += t1;
  
          num2 = sparse_binary_points[i2].id.size();
          for (p2=0; p2<num2; p2++)
              w[sparse_binary_points[i2].id[p2]] += t2;
      }
      else if (is_sparse_data && !is_binary) {
          int p1,num1,p2,num2;
  
          num1 = sparse_points[i1].id.size();
          for (p1=0; p1<num1; p1++)
              w[sparse_points[i1].id[p1]] += 
                  t1 * sparse_points[i1].val[p1];
  
          num2 = sparse_points[i2].id.size();
          for (p2=0; p2<num2; p2++)
              w[sparse_points[i2].id[p2]] +=
                  t2 * sparse_points[i2].val[p2];
      }
      else
          for (int i=0; i<d; i++)
              w[i] += dense_points[i1][i] * t1 + dense_points[i2][i] * t2;
  }
#line 1102 "smo.w"

  
#line 1262 "smo.w"
  {
      float t1 = y1 * (a1-alph1);
      float t2 = y2 * (a2-alph2);
  
      for (int i=0; i<end_support_i; i++)
      if (0 < alph[i] && alph[i] < C)
          error_cache[i] +=  t1 * kernel_func(i1,i) + t2 * kernel_func(i2,i)
                            - delta_b;
      error_cache[i1] = 0.;
      error_cache[i2] = 0.;
  }
#line 1103 "smo.w"


  alph[i1] = a1;  /* Store a1 in the alpha array.*/
  alph[i2] = a2;  /* Store a2 in the alpha array.*/

  return 1;
}
#line 1286 "smo.w"
float learned_func_linear_sparse_binary(int k) {
  float s = 0.;

  for (int i=0; i<sparse_binary_points[k].id.size(); i++)
    s += w[sparse_binary_points[k].id[i]];

  s -= b;
  return s;
}
#line 1297 "smo.w"
float learned_func_linear_sparse_nonbinary(int k) {
  float s = 0.;

  for (int i=0; i<sparse_points[k].id.size(); i++)
  {
      int j = sparse_points[k].id[i];
      float v = sparse_points[k].val[i];
      s += w[j] * v;
  }
  s -= b;
  return s;
}
#line 1311 "smo.w"
float learned_func_linear_dense(int k) {
  float s = 0.;

  for (int i=0; i<d; i++)
    s += w[i] * dense_points[k][i];

  s -= b;
  return s;
}
#line 1322 "smo.w"
float learned_func_nonlinear(int k) {
  float s = 0.;
  for (int i=0; i<end_support_i; i++)
      if (alph[i] > 0)
          s += alph[i]*target[i]*kernel_func(i,k);
  s -= b;
  return s;
}
#line 1366 "smo.w"

float dot_product_sparse_binary(int i1, int i2)
{
  int p1=0, p2=0, dot=0;
  int num1 = sparse_binary_points[i1].id.size();
  int num2 = sparse_binary_points[i2].id.size();

  while (p1 < num1 && p2 < num2) {
      int a1 = sparse_binary_points[i1].id[p1];
      int a2 = sparse_binary_points[i2].id[p2];
      if (a1 == a2) {
          dot++;
          p1++;
          p2++;
        }
      else if (a1 > a2)
        p2++;
      else
        p1++;
    }
  return (float)dot;
}
#line 1390 "smo.w"

float dot_product_sparse_nonbinary(int i1, int i2)
{
  int p1=0, p2=0;
  float dot = 0.;
  int num1 = sparse_points[i1].id.size();
  int num2 = sparse_points[i2].id.size();

  while (p1 < num1 && p2 < num2) {