svm.cs

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

					{
						if (- G[i] > Gmax3)
						{
							Gmax3 = - G[i];
							Gmax3_idx = i;
						}
					}
					if (!is_lower_bound(i))
					// d = -1
					{
						if (G[i] > Gmax4)
						{
							Gmax4 = G[i];
							Gmax4_idx = i;
						}
					}
				}
			}
			
			if (System.Math.Max(Gmax1 + Gmax2, Gmax3 + Gmax4) < eps)
				return 1;
			
			if (Gmax1 + Gmax2 > Gmax3 + Gmax4)
			{
				working_set[0] = Gmax1_idx;
				working_set[1] = Gmax2_idx;
			}
			else
			{
				working_set[0] = Gmax3_idx;
				working_set[1] = Gmax4_idx;
			}
			return 0;
		}
		
		internal override void  do_shrinking()
		{
			double Gmax1 = - INF; // max { -grad(f)_i * d | y_i = +1, d = +1 }
			double Gmax2 = - INF; // max { -grad(f)_i * d | y_i = +1, d = -1 }
			double Gmax3 = - INF; // max { -grad(f)_i * d | y_i = -1, d = +1 }
			double Gmax4 = - INF; // max { -grad(f)_i * d | y_i = -1, d = -1 }
			
			int k;
			for (k = 0; k < active_size; k++)
			{
				if (!is_upper_bound(k))
				{
					if (y[k] == + 1)
					{
						if (- G[k] > Gmax1)
							Gmax1 = - G[k];
					}
					else if (- G[k] > Gmax3)
						Gmax3 = - G[k];
				}
				if (!is_lower_bound(k))
				{
					if (y[k] == + 1)
					{
						if (G[k] > Gmax2)
							Gmax2 = G[k];
					}
					else if (G[k] > Gmax4)
						Gmax4 = G[k];
				}
			}
			
			double Gm1 = - Gmax2;
			double Gm2 = - Gmax1;
			double Gm3 = - Gmax4;
			double Gm4 = - Gmax3;
			
			for (k = 0; k < active_size; k++)
			{
				if (is_lower_bound(k))
				{
					if (y[k] == + 1)
					{
						if (- G[k] >= Gm1)
							continue;
					}
					else if (- G[k] >= Gm3)
						continue;
				}
				else if (is_upper_bound(k))
				{
					if (y[k] == + 1)
					{
						if (G[k] >= Gm2)
							continue;
					}
					else if (G[k] >= Gm4)
						continue;
				}
				else
					continue;
				
				--active_size;
				swap_index(k, active_size);
				--k; // look at the newcomer
			}
			
			// unshrink, check all variables again before final iterations
			
			if (unshrinked || System.Math.Max(- (Gm1 + Gm2), - (Gm3 + Gm4)) > eps * 10)
				return ;
			
			unshrinked = true;
			reconstruct_gradient();
			
			for (k = l - 1; k >= active_size; k--)
			{
				if (is_lower_bound(k))
				{
					if (y[k] == + 1)
					{
						if (- G[k] < Gm1)
							continue;
					}
					else if (- G[k] < Gm3)
						continue;
				}
				else if (is_upper_bound(k))
				{
					if (y[k] == + 1)
					{
						if (G[k] < Gm2)
							continue;
					}
					else if (G[k] < Gm4)
						continue;
				}
				else
					continue;
				
				swap_index(k, active_size);
				active_size++;
				++k; // look at the newcomer
			}
		}
		
		internal override double calculate_rho()
		{
			int nr_free1 = 0, nr_free2 = 0;
			double ub1 = INF, ub2 = INF;
			double lb1 = - INF, lb2 = - INF;
			double sum_free1 = 0, sum_free2 = 0;
			
			for (int i = 0; i < active_size; i++)
			{
				if (y[i] == + 1)
				{
					if (is_lower_bound(i))
						ub1 = System.Math.Min(ub1, G[i]);
					else if (is_upper_bound(i))
						lb1 = System.Math.Max(lb1, G[i]);
					else
					{
						++nr_free1;
						sum_free1 += G[i];
					}
				}
				else
				{
					if (is_lower_bound(i))
						ub2 = System.Math.Min(ub2, G[i]);
					else if (is_upper_bound(i))
						lb2 = System.Math.Max(lb2, G[i]);
					else
					{
						++nr_free2;
						sum_free2 += G[i];
					}
				}
			}
			
			double r1, r2;
			if (nr_free1 > 0)
				r1 = sum_free1 / nr_free1;
			else
				r1 = (ub1 + lb1) / 2;
			
			if (nr_free2 > 0)
				r2 = sum_free2 / nr_free2;
			else
				r2 = (ub2 + lb2) / 2;
			
			si.r = (r1 + r2) / 2;
			return (r1 - r2) / 2;
		}
	}
	
	//
	// Q matrices for various formulations
	//
	class SVC_Q:Kernel
	{
		//UPGRADE_NOTE: Final was removed from the declaration of 'y '. 'ms-help://MS.VSCC.2003/commoner/redir/redirect.htm?keyword="jlca1003_3"'
		private sbyte[] y;
		//UPGRADE_NOTE: Final was removed from the declaration of 'cache '. 'ms-help://MS.VSCC.2003/commoner/redir/redirect.htm?keyword="jlca1003_3"'
		private Cache cache;
		
		internal SVC_Q(svm_problem prob, svm_parameter param, sbyte[] y_):base(prob.l, prob.x, param)
		{
			y = new sbyte[y_.Length];
			y_.CopyTo(y, 0);
			//UPGRADE_WARNING: Data types in Visual C# might be different.  Verify the accuracy of narrowing conversions. 'ms-help://MS.VSCC.2003/commoner/redir/redirect.htm?keyword="jlca1042_3"'
			cache = new Cache(prob.l, (int) (param.cache_size * (1 << 20)));
		}
		
		internal override float[] get_Q(int i, int len)
		{
			float[][] data = new float[1][];
			int start;
			if ((start = cache.get_data(i, data, len)) < len)
			{
				for (int j = start; j < len; j++)
				{
					//UPGRADE_WARNING: Data types in Visual C# might be different.  Verify the accuracy of narrowing conversions. 'ms-help://MS.VSCC.2003/commoner/redir/redirect.htm?keyword="jlca1042_3"'
					data[0][j] = (float) (y[i] * y[j] * kernel_function(i, j));
				}
			}
			return data[0];
		}
		
		internal override void  swap_index(int i, int j)
		{
			cache.swap_index(i, j);
			base.swap_index(i, j);
			do 
			{
				sbyte _ = y[i]; y[i] = y[j]; y[j] = _;
			}
			while (false);
		}
	}
	
	class ONE_CLASS_Q:Kernel
	{
		//UPGRADE_NOTE: Final was removed from the declaration of 'cache '. 'ms-help://MS.VSCC.2003/commoner/redir/redirect.htm?keyword="jlca1003_3"'
		private Cache cache;
		
		internal ONE_CLASS_Q(svm_problem prob, svm_parameter param):base(prob.l, prob.x, param)
		{
			//UPGRADE_WARNING: Data types in Visual C# might be different.  Verify the accuracy of narrowing conversions. 'ms-help://MS.VSCC.2003/commoner/redir/redirect.htm?keyword="jlca1042_3"'
			cache = new Cache(prob.l, (int) (param.cache_size * (1 << 20)));
		}
		
		internal override float[] get_Q(int i, int len)
		{
			float[][] data = new float[1][];
			int start;
			if ((start = cache.get_data(i, data, len)) < len)
			{
				for (int j = start; j < len; j++)
				{
					//UPGRADE_WARNING: Data types in Visual C# might be different.  Verify the accuracy of narrowing conversions. 'ms-help://MS.VSCC.2003/commoner/redir/redirect.htm?keyword="jlca1042_3"'
					data[0][j] = (float) kernel_function(i, j);
				}
			}
			return data[0];
		}
		
		internal override void  swap_index(int i, int j)
		{
			cache.swap_index(i, j);
			base.swap_index(i, j);
		}
	}
	
	class SVR_Q:Kernel
	{
		//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;
		//UPGRADE_NOTE: Final was removed from the declaration of 'cache '. 'ms-help://MS.VSCC.2003/commoner/redir/redirect.htm?keyword="jlca1003_3"'
		private Cache cache;
		//UPGRADE_NOTE: Final was removed from the declaration of 'sign '. 'ms-help://MS.VSCC.2003/commoner/redir/redirect.htm?keyword="jlca1003_3"'
		private sbyte[] sign;
		//UPGRADE_NOTE: Final was removed from the declaration of 'index '. 'ms-help://MS.VSCC.2003/commoner/redir/redirect.htm?keyword="jlca1003_3"'
		private int[] index;
		private int next_buffer;
		private float[][] buffer;
		
		internal SVR_Q(svm_problem prob, svm_parameter param):base(prob.l, prob.x, param)
		{
			l = prob.l;
			//UPGRADE_WARNING: Data types in Visual C# might be different.  Verify the accuracy of narrowing conversions. 'ms-help://MS.VSCC.2003/commoner/redir/redirect.htm?keyword="jlca1042_3"'
			cache = new Cache(l, (int) (param.cache_size * (1 << 20)));
			sign = new sbyte[2 * l];
			index = new int[2 * l];
			for (int k = 0; k < l; k++)
			{
				sign[k] = 1;
				sign[k + l] = - 1;
				index[k] = k;
				index[k + l] = k;
			}
			buffer = new float[2][];
			for (int i = 0; i < 2; i++)
			{
				buffer[i] = new float[2 * l];
			}
			next_buffer = 0;
		}
		
		internal override void  swap_index(int i, int j)
		{
			do 
			{
				sbyte _ = sign[i]; sign[i] = sign[j]; sign[j] = _;
			}
			while (false);
			do 
			{
				int _ = index[i]; index[i] = index[j]; index[j] = _;
			}
			while (false);
		}
		
		internal override float[] get_Q(int i, int len)
		{
			float[][] data = new float[1][];
			int real_i = index[i];
			if (cache.get_data(real_i, data, l) < l)
			{
				for (int j = 0; j < l; j++)
				{
					//UPGRADE_WARNING: Data types in Visual C# might be different.  Verify the accuracy of narrowing conversions. 'ms-help://MS.VSCC.2003/commoner/redir/redirect.htm?keyword="jlca1042_3"'
					data[0][j] = (float) kernel_function(real_i, j);
				}
			}
			
			// reorder and copy
			float[] buf = buffer[next_buffer];
			next_buffer = 1 - next_buffer;
			sbyte si = sign[i];
			for (int j = 0; j < len; j++)
				buf[j] = si * sign[j] * data[0][index[j]];
			return buf;
		}
	}
	
	public class svm
	{
		//
		// construct and solve various formulations
		//
		private static void  solve_c_svc(svm_problem prob, svm_parameter param, double[] alpha, Solver.SolutionInfo si, double Cp, double Cn)
		{
			int l = prob.l;
			double[] minus_ones = new double[l];
			sbyte[] y = new sbyte[l];
			
			int i;
			
			for (i = 0; i < l; i++)
			{
				alpha[i] = 0;
				minus_ones[i] = - 1;
				if (prob.y[i] > 0)
					y[i] = (sbyte) (+ 1);
				else
					y[i] = - 1;
			}
			
			Solver s = new Solver();
			s.Solve(l, new SVC_Q(prob, param, y), minus_ones, y, alpha, Cp, Cn, param.eps, si, param.shrinking);
			
			double sum_alpha = 0;
			for (i = 0; i < l; i++)
				sum_alpha += alpha[i];
			
			if (Cp == Cn)
				System.Console.Out.Write("nu = " + sum_alpha / (Cp * prob.l) + "\n");
			
			for (i = 0; i < l; i++)
				alpha[i] *= y[i];
		}
		
		private static void  solve_nu_svc(svm_problem prob, svm_parameter param, double[] alpha, Solver.SolutionInfo si)
		{
			int i;
			int l = prob.l;
			double nu = param.nu;
			
			sbyte[] y = new sbyte[l];
			
			for (i = 0; i < l; i++)
				if (prob.y[i] > 0)
					y[i] = (sbyte) (+ 1);
				else
					y[i] = - 1;
			
			double sum_pos = nu * l / 2;
			double sum_neg = nu * l / 2;
			
			for (i = 0; i < l; i++)
				if (y[i] == + 1)
				{
					alpha[i] = System.Math.Min(1.0, sum_pos);
					sum_pos -= alpha[i];
				}
				else
				{
					alpha[i] = System.Math.Min(1.0, sum_neg);
					sum_neg -= alpha[i];
				}
			
			double[] zeros = new double[l];
			
			for (i = 0; i < l; i++)
				zeros[i] = 0;
			
			Solver_NU s = new Solver_NU();
			s.Solve(l, new SVC_Q(prob, param, y), zeros, y, alpha, 1.0, 1.0, param.eps, si, param.shrinking);
			double r = si.r;
			
			System.Console.Out.Write("C = " + 1 / r + "\n");
			
			for (i = 0; i < l; i++)
				alpha[i] *= y[i] / r;
			
			si.rho /= r;
			si.obj /= (r * r);
			si.upper_bound_p = 1 / r;
			si.upper_bound_n = 1 / r;
		}
		
		private static void  solve_one_class(svm_problem prob, svm_parameter param, double[] alpha, Solver.SolutionInfo si)
		{
			int l = prob.l;
			double[] zeros = new double[l];
			sbyte[] ones = new sbyte[l];
			int i;
			
			//UPGRADE_WARNING: Data types in Visual C# might be different.  Verify the accuracy of narrowing conversions. 'ms-help://MS.VSCC.2003/commoner/redir/redirect.htm?keyword="jlca1042_3"'
			int n = (int) (param.nu * prob.l); // # of alpha's at upper bound
			
			for (i = 0; i < n; i++)
				alpha[i] = 1;
			alpha[n] = param.nu * prob.l - n;
			for (i = n + 1; i < l; i++)
				alpha[i] = 0;
			
			for (i = 0; i < l; i++)
			{
				zeros[i] = 0;
				ones[i] = 1;
			}
			
			Solver s = new Solver();
			s.Solve(l, new ONE_CLASS_Q(prob, param), zeros, ones, alpha, 1.0, 1.0, param.eps, si, param.shrinking);
		}
		
		private static void  solve_epsilon_svr(svm_problem prob, svm_parameter param, double[] alpha, Solver.SolutionInfo si)
		{
			int l = prob.l;
			double[] alpha2 = new double[2 * l];
			double[] linear_term = new double[2 * l];
			sbyte[] y = new sbyte[2 * l];
			int i;
			
			for (i = 0; i < l; i++)
			{
				alpha2[i] = 0;
				linear_term[i] = param.p - prob.y[i];
				y[i] = 1;
				
				alpha2[i + l] = 0;
				linear_term[i + l] = param.p + prob.y[i];
				y[i + l] = - 1;
			}
			
			Solver s = new Solver();
			s.Solve(2 * l, new SVR_Q(prob, param), linear_term, y, alpha2, param.C, param.C, param.eps, si, param.shrinking);
			
			double sum_alpha = 0;
			for (i = 0; i < l; i++)
			{
				alpha[i] = alpha2[i] - alpha2[i + l];
				sum_alpha += System.Math.Abs(alpha[i]);
			}
			System.Console.Out.Write("nu = " + sum_alpha / (param.C * l) + "\n");
		}
		
		private static void  solve_nu_svr(svm_problem prob, svm_parameter param, double[] alpha, Solver.SolutionInfo si)
		{
			int l = prob.l;
			double C = param.C;
			double[] alpha2 = new double[2 * l];
			double[] linear_term = new double[2 * l];
			sbyte[] y = new sbyte[2 * l];
			int i;