Package com.nr.ci

Source Code of com.nr.ci.Svm

package com.nr.ci;
import static java.lang.Math.*;
import static com.nr.NRUtil.*;
import com.nr.ran.Ran;
import com.nr.sort.Indexx;

/**
* Support Vector Machines
* Copyright (C) Numerical Recipes Software 1986-2007
* Java translation Copyright (C) Huang Wen Hui 2012
*
* @author hwh
*
*/

// class for solving SVM problems by the SOR method
public class Svm {
  private Svmgenkernel gker;  // Reference bound to user's kernel (and data)
  private int m, fnz, fub, niter;
  private double[] alph, alphold;   // Vectors of a's before and after a step
  private Ran ran;   // Random number generator
  private boolean alphinit;
  private double dalph;   // Change in norm of the a's in one step
 
  // constructor binds the user's kernel and allocates storage
  public Svm(final Svmgenkernel inker){
    gker = inker;
    m = gker.y.length;
    alph = new double[m];
    alphold= new double[m];
    ran = new Ran(21);
    alphinit = false;
  }
 
  // Perform one group of relaxation steps: a single step over all the a's, and multiple steps over only the interior a's
  public double relax(final double lambda, final double om) {
    int iter,j,jj,k,kk;
    double sum;   // index when a's are sorted by value
    double[] pinsum = new double[m]// stored sums over noninterior variables
    if (alphinit == false) {   // start all a's at 0
      for (j=0; j<m; j++) alph[j] = 0.;
      alphinit = true;
    }
    alphold = alph;   // save old a's
    // here begins the relaxation pass over all the a's
    Indexx x = new Indexx(alph)// sort a's, then find first nonzero one
    for (fnz=0; fnz<m; fnz++) if (alph[x.indx[fnz]] != 0.) break;
    for (j=fnz; j<m-2; j++) {  // randomly permute all the nonzero a's
      k = j + (ran.int32p() % (m-j));
      swap(x.indx,j,k);
    }
    for (jj=0; jj<m; jj++) {   // Main loop over a's
      j = x.indx[jj];
      sum = 0.;
      for (kk=fnz; kk<m; kk++) {  // Sums start with first nonzero
        k = x.indx[kk];
        sum += (gker.ker[j][k] + 1.)*gker.y[k]*alph[k];
      }
      alph[j] = alph[j] - (om/(gker.ker[j][j]+1.))*(gker.y[j]*sum-1.);
      alph[j] = max(0.,min(lambda,alph[j]));   // Projection operator
      if (jj < fnz && alph[j]!=0) {
        --fnz;
        //SWAP(x.indx[--fnz],x.indx[jj]);
        swap(x.indx, fnz, jj);
      // (Above) MAke an \alpha active if it becomes nonzero
    }
   
    // Here begins the relaxation passes over the interior \alpha's
    Indexx y = new Indexx(alph);   // Sort. Identify interior \alpha's
    for (fnz=0; fnz<m; fnz++) if (alph[y.indx[fnz]] != 0.) break;
    for (fub=fnz; fub<m; fub++) if (alph[y.indx[fub]] == lambda) break;
    for (j=fnz; j<fub-2; j++) {  // Permute
      k = j + (ran.int32p() % (fub-j));
      swap(y.indx,j,k);
    }
   
    for (jj=fnz; jj<fub; jj++) {  // Compute sums over pinned \alpha's just once
      j = y.indx[jj];
      sum = 0.;
      for (kk=fub; kk<m; kk++) {
        k = y.indx[kk];
        sum += (gker.ker[j][k] + 1.)*gker.y[k]*alph[k];
      }
      pinsum[jj] = sum;
    }
    niter = max((int)(0.5*(m+1.0)*(m-fnz+1.0)/(SQR(fub-fnz+1.0))),1);
   
    // Calculate a numer of iterations that will take about half as long as the full pass just completed
    for (iter=0; iter<niter; iter++) {  // Main loop over \alpha's
      for (jj=fnz; jj<fub; jj++) {
        j = y.indx[jj];
        sum = pinsum[jj];
        for (kk=fnz; kk<fub; kk++) {
          k = y.indx[kk];
          sum += (gker.ker[j][k] + 1.)*gker.y[k]*alph[k];
        }
        alph[j] = alph[j] - (om/(gker.ker[j][j]+1.))*(gker.y[j]*sum-1.);
        alph[j] = max(0.,min(lambda,alph[j]));
      }  
    }
    dalph = 0.// Return change in norm of \alpha vector
    for (j=0;j<m;j++) dalph += SQR(alph[j]-alphold[j]);
    return sqrt(dalph);
  }
 
 
  // Call only after convergence via repeated calls to relax. Returns the decision rule f(x) for data point k
  public double predict(final int k) {
    double sum = 0.;
    for (int j=0; j<m; j++) sum += alph[j]*gker.y[j]*(gker.ker[j][k]+1.0);
    return sum;
  }
 
  // Call only after convergence via repeated calls to relax. Returns the decision rule f(x) for an arbitrary feature vector
  public double predict(final double[] x) {
    double sum = 0.;
    for (int j=0; j<m; j++) sum += alph[j]*gker.y[j]*(gker.kernel(j,x)+1.0);
    return sum;
  }
}
TOP

Related Classes of com.nr.ci.Svm

TOP
Copyright © 2018 www.massapi.com. All rights reserved.
All source code are property of their respective owners. Java is a trademark of Sun Microsystems, Inc and owned by ORACLE Inc. Contact coftware#gmail.com.