svm.cs

这是C#版本开发的SVM类库包,适合不同爱好的同学学习.

/*
 * Conversion notes:
 * Using JLCA 3.0, the only problem was with the save and loads. I changed them both to be
 * StreamR/W around a FileStream. Originally, the save was a BinaryWriter over a FileStream
 * and the Reader was a StreamReader over another StreamReader.
 * In the Java code, it's a DataOutputStream around a FileOutputStream and a BufferedReader
 * around a FileReader.
 */

using System;
namespace libsvm
{
	
	//
	// Kernel Cache
	//
	// l is the number of total data items
	// size is the cache size limit in bytes
	//
	class Cache
	{
		//UPGRADE_NOTE: Final was removed from the declaration of 'l '. 'ms-help://MS.VSCC.2003/commoner/redir/redirect.htm?keyword="jlca1003_3"'
		private int l;
		private int size;
		//UPGRADE_NOTE: Field 'EnclosingInstance' was added to class 'head_t' to access its enclosing instance. 'ms-help://MS.VSCC.2003/commoner/redir/redirect.htm?keyword="jlca1019_3"'
		private sealed class head_t
		{
			public head_t(Cache enclosingInstance)
			{
				InitBlock(enclosingInstance);
			}
			private void  InitBlock(Cache enclosingInstance)
			{
				this.enclosingInstance = enclosingInstance;
			}
			private Cache enclosingInstance;
			public Cache Enclosing_Instance
			{
				get
				{
					return enclosingInstance;
				}
				
			}
			internal head_t prev, next; // a cicular list
			internal float[] data;
			internal int len; // data[0,len) is cached in this entry
		}
		//UPGRADE_NOTE: Final was removed from the declaration of 'head '. 'ms-help://MS.VSCC.2003/commoner/redir/redirect.htm?keyword="jlca1003_3"'
		private head_t[] head;
		private head_t lru_head;
		
		internal Cache(int l_, int size_)
		{
			l = l_;
			size = size_;
			head = new head_t[l];
			for (int i = 0; i < l; i++)
				head[i] = new head_t(this);
			size /= 4;
			size -= l * (16 / 4); // sizeof(head_t) == 16
			lru_head = new head_t(this);
			lru_head.next = lru_head.prev = lru_head;
		}
		
		private void  lru_delete(head_t h)
		{
			// delete from current location
			h.prev.next = h.next;
			h.next.prev = h.prev;
		}
		
		private void  lru_insert(head_t h)
		{
			// insert to last position
			h.next = lru_head;
			h.prev = lru_head.prev;
			h.prev.next = h;
			h.next.prev = h;
		}
		
		// request data [0,len)
		// return some position p where [p,len) need to be filled
		// (p >= len if nothing needs to be filled)
		// java: simulate pointer using single-element array
		internal virtual int get_data(int index, float[][] data, int len)
		{
			head_t h = head[index];
			if (h.len > 0)
				lru_delete(h);
			int more = len - h.len;
			
			if (more > 0)
			{
				// free old space
				while (size < more)
				{
					head_t old = lru_head.next;
					lru_delete(old);
					size += old.len;
					old.data = null;
					old.len = 0;
				}
				
				// allocate new space
				float[] new_data = new float[len];
				if (h.data != null)
					Array.Copy(h.data, 0, new_data, 0, h.len);
				h.data = new_data;
				size -= more;
				do 
				{
					int _ = h.len; h.len = len; len = _;
				}
				while (false);
			}
			
			lru_insert(h);
			data[0] = h.data;
			return len;
		}
		
		internal virtual void  swap_index(int i, int j)
		{
			if (i == j)
				return ;
			
			if (head[i].len > 0)
				lru_delete(head[i]);
			if (head[j].len > 0)
				lru_delete(head[j]);
			do 
			{
				float[] _ = head[i].data; head[i].data = head[j].data; head[j].data = _;
			}
			while (false);
			do 
			{
				int _ = head[i].len; head[i].len = head[j].len; head[j].len = _;
			}
			while (false);
			if (head[i].len > 0)
				lru_insert(head[i]);
			if (head[j].len > 0)
				lru_insert(head[j]);
			
			if (i > j)
				do 
				{
					int _ = i; i = j; j = _;
				}
				while (false);
			for (head_t h = lru_head.next; h != lru_head; h = h.next)
			{
				if (h.len > i)
				{
					if (h.len > j)
						do 
						{
							float _ = h.data[i]; h.data[i] = h.data[j]; h.data[j] = _;
						}
						while (false);
					else
					{
						// give up
						lru_delete(h);
						size += h.len;
						h.data = null;
						h.len = 0;
					}
				}
			}
		}
	}
	
	//
	// Kernel evaluation
	//
	// the static method k_function is for doing single kernel evaluation
	// the constructor of Kernel prepares to calculate the l*l kernel matrix
	// the member function get_Q is for getting one column from the Q Matrix
	//
	abstract class Kernel
	{
		private svm_node[][] x;
		//UPGRADE_NOTE: Final was removed from the declaration of 'x_square '. 'ms-help://MS.VSCC.2003/commoner/redir/redirect.htm?keyword="jlca1003_3"'
		private double[] x_square;
		
		// svm_parameter
		//UPGRADE_NOTE: Final was removed from the declaration of 'kernel_type '. 'ms-help://MS.VSCC.2003/commoner/redir/redirect.htm?keyword="jlca1003_3"'
		private int kernel_type;
		//UPGRADE_NOTE: Final was removed from the declaration of 'degree '. 'ms-help://MS.VSCC.2003/commoner/redir/redirect.htm?keyword="jlca1003_3"'
		private double degree;
		//UPGRADE_NOTE: Final was removed from the declaration of 'gamma '. 'ms-help://MS.VSCC.2003/commoner/redir/redirect.htm?keyword="jlca1003_3"'
		private double gamma;
		//UPGRADE_NOTE: Final was removed from the declaration of 'coef0 '. 'ms-help://MS.VSCC.2003/commoner/redir/redirect.htm?keyword="jlca1003_3"'
		private double coef0;
		
		internal abstract float[] get_Q(int column, int len);
		
		internal virtual void  swap_index(int i, int j)
		{
			do 
			{
				svm_node[] _ = x[i]; x[i] = x[j]; x[j] = _;
			}
			while (false);
			if (x_square != null)
				do 
				{
					double _ = x_square[i]; x_square[i] = x_square[j]; x_square[j] = _;
				}
				while (false);
		}
		
		private static double tanh(double x)
		{
			double e = System.Math.Exp(x);
			return 1.0 - 2.0 / (e * e + 1);
		}
		
		internal virtual double kernel_function(int i, int j)
		{
			switch (kernel_type)
			{
				
				case svm_parameter.LINEAR: 
					return dot(x[i], x[j]);
				
				case svm_parameter.POLY: 
					return System.Math.Pow(gamma * dot(x[i], x[j]) + coef0, degree);
				
				case svm_parameter.RBF: 
					return System.Math.Exp((- gamma) * (x_square[i] + x_square[j] - 2 * dot(x[i], x[j])));
				
				case svm_parameter.SIGMOID: 
					return tanh(gamma * dot(x[i], x[j]) + coef0);
				
				default: 
					return 0; // java
				
			}
		}
		
		internal Kernel(int l, svm_node[][] x_, svm_parameter param)
		{
			this.kernel_type = param.kernel_type;
			this.degree = param.degree;
			this.gamma = param.gamma;
			this.coef0 = param.coef0;
			
			x = (svm_node[][]) x_.Clone();
			
			if (kernel_type == svm_parameter.RBF)
			{
				x_square = new double[l];
				for (int i = 0; i < l; i++)
					x_square[i] = dot(x[i], x[i]);
			}
			else
				x_square = null;
		}
		
		internal static double dot(svm_node[] x, svm_node[] y)
		{
			double sum = 0;
			int xlen = x.Length;
			int ylen = y.Length;
			int i = 0;
			int j = 0;
			while (i < xlen && j < ylen)
			{
				if (x[i].index == y[j].index)
					sum += x[i++].value_Renamed * y[j++].value_Renamed;
				else
				{
					if (x[i].index > y[j].index)
						++j;
					else
						++i;
				}
			}
			return sum;
		}
		
		internal static double k_function(svm_node[] x, svm_node[] y, svm_parameter param)
		{
			switch (param.kernel_type)
			{
				
				case svm_parameter.LINEAR: 
					return dot(x, y);
				
				case svm_parameter.POLY: 
					return System.Math.Pow(param.gamma * dot(x, y) + param.coef0, param.degree);
				
				case svm_parameter.RBF: 
				{
					double sum = 0;
					int xlen = x.Length;
					int ylen = y.Length;
					int i = 0;
					int j = 0;
					while (i < xlen && j < ylen)
					{
						if (x[i].index == y[j].index)
						{
							double d = x[i++].value_Renamed - y[j++].value_Renamed;
							sum += d * d;
						}
						else if (x[i].index > y[j].index)
						{
							sum += y[j].value_Renamed * y[j].value_Renamed;
							++j;
						}
						else
						{
							sum += x[i].value_Renamed * x[i].value_Renamed;
							++i;
						}
					}
					
					while (i < xlen)
					{
						sum += x[i].value_Renamed * x[i].value_Renamed;
						++i;
					}
					
					while (j < ylen)
					{
						sum += y[j].value_Renamed * y[j].value_Renamed;
						++j;
					}
					
					return System.Math.Exp((- param.gamma) * sum);
				}
				
				case svm_parameter.SIGMOID: 
					return tanh(param.gamma * dot(x, y) + param.coef0);
				
				default: 
					return 0; // java
				
			}
		}
	}
	
	// Generalized SMO+SVMlight algorithm
	// Solves:
	//
	//	min 0.5(\alpha^T Q \alpha) + b^T \alpha
	//
	//		y^T \alpha = \delta
	//		y_i = +1 or -1
	//		0 <= alpha_i <= Cp for y_i = 1
	//		0 <= alpha_i <= Cn for y_i = -1
	//
	// Given:
	//
	//	Q, b, y, Cp, Cn, and an initial feasible point \alpha
	//	l is the size of vectors and matrices
	//	eps is the stopping criterion
	//
	// solution will be put in \alpha, objective value will be put in obj
	//
	class Solver
	{
		internal int active_size;
		internal sbyte[] y;
		internal double[] G; // gradient of objective function
		internal const sbyte LOWER_BOUND = 0;
		internal const sbyte UPPER_BOUND = 1;
		internal const sbyte FREE = 2;
		internal sbyte[] alpha_status; // LOWER_BOUND, UPPER_BOUND, FREE
		internal double[] alpha;
		internal Kernel Q;
		internal double eps;
		internal double Cp, Cn;
		internal double[] b;
		internal int[] active_set;
		internal double[] G_bar; // gradient, if we treat free variables as 0
		internal int l;
		internal bool unshrinked; // XXX
		
		//UPGRADE_NOTE: Final was removed from the declaration of 'INF '. 'ms-help://MS.VSCC.2003/commoner/redir/redirect.htm?keyword="jlca1003_3"'
		internal static readonly double INF = System.Double.PositiveInfinity;
		
		internal virtual double get_C(int i)
		{
			return (y[i] > 0)?Cp:Cn;
		}
		internal virtual void  update_alpha_status(int i)
		{
			if (alpha[i] >= get_C(i))
				alpha_status[i] = UPPER_BOUND;
			else if (alpha[i] <= 0)
				alpha_status[i] = LOWER_BOUND;
			else
				alpha_status[i] = FREE;
		}
		internal virtual bool is_upper_bound(int i)
		{
			return alpha_status[i] == UPPER_BOUND;
		}
		internal virtual bool is_lower_bound(int i)
		{
			return alpha_status[i] == LOWER_BOUND;
		}
		internal virtual bool is_free(int i)
		{
			return alpha_status[i] == FREE;
		}
		
		// java: information about solution except alpha,
		// because we cannot return multiple values otherwise...
		internal class SolutionInfo
		{
			internal double obj;
			internal double rho;
			internal double upper_bound_p;
			internal double upper_bound_n;
			internal double r; // for Solver_NU
		}
		
		internal virtual void  swap_index(int i, int j)
		{
			Q.swap_index(i, j);
			do 
			{
				sbyte _ = y[i]; y[i] = y[j]; y[j] = _;
			}
			while (false);
			do 
			{
				double _ = G[i]; G[i] = G[j]; G[j] = _;
			}
			while (false);
			do 
			{
				sbyte _ = alpha_status[i]; alpha_status[i] = alpha_status[j]; alpha_status[j] = _;
			}
			while (false);
			do 
			{
				double _ = alpha[i]; alpha[i] = alpha[j]; alpha[j] = _;
			}
			while (false);
			do 
			{
				double _ = b[i]; b[i] = b[j]; b[j] = _;
			}
			while (false);
			do 
			{
				int _ = active_set[i]; active_set[i] = active_set[j]; active_set[j] = _;
			}
			while (false);
			do 
			{
				double _ = G_bar[i]; G_bar[i] = G_bar[j]; G_bar[j] = _;
			}
			while (false);
		}
		
		internal virtual void  reconstruct_gradient()
		{
			// reconstruct inactive elements of G from G_bar and free variables
			
			if (active_size == l)
				return ;
			
			int i;
			for (i = active_size; i < l; i++)
				G[i] = G_bar[i] + b[i];
			
			for (i = 0; i < active_size; i++)
				if (is_free(i))
				{
					float[] Q_i = Q.get_Q(i, l);
					double alpha_i = alpha[i];
					for (int j = active_size; j < l; j++)
						G[j] += alpha_i * Q_i[j];
				}
		}
		
		internal virtual void  Solve(int l, Kernel Q, double[] b_, sbyte[] y_, double[] alpha_, double Cp, double Cn, double eps, SolutionInfo si, int shrinking)
		{
			this.l = l;
			this.Q = Q;
			b = new double[b_.Length];
			b_.CopyTo(b, 0);
			y = new sbyte[y_.Length];
			y_.CopyTo(y, 0);