if (grad_diff > 0)					{						double obj_diff; 						double quad_coef=Q_i[i]+QD[j]+2*y[i]*Q_i[j];						if (quad_coef > 0)							obj_diff = -(grad_diff*grad_diff)/quad_coef;						else							obj_diff = -(grad_diff*grad_diff)/TAU;							if (obj_diff <= obj_diff_min)						{							Gmin_idx=j;							obj_diff_min = obj_diff;						}					}				}			}		}		if(Gmax+Gmax2 < eps)			return 1;		working_set[0] = Gmax_idx;		working_set[1] = Gmin_idx;		return 0;	}	// return 1 if already optimal, return 0 otherwise	int max_violating_pair(int[] working_set)	{		// return i,j which maximize -grad(f)^T d , under constraint		// if alpha_i == C, d != +1		// if alpha_i == 0, d != -1		double Gmax1 = -INF;		// max { -y_i * grad(f)_i | i in I_up(\alpha) }		int Gmax1_idx = -1;		int Gmax2_idx = -1;		double Gmax2 = -INF;		// max { y_i * grad(f)_i | i in I_low(\alpha) }		for(int i=0;i<active_size;i++)		{			if(y[i]==+1)	// y = +1			{				if(!is_upper_bound(i))	// d = +1				{					if(-G[i] >= Gmax1)					{						Gmax1 = -G[i];						Gmax1_idx = i;					}				}				if(!is_lower_bound(i))	// d = -1				{					if(G[i] >= Gmax2)					{						Gmax2 = G[i];						Gmax2_idx = i;					}				}			}			else		// y = -1			{				if(!is_upper_bound(i))	// d = +1				{					if(-G[i] >= Gmax2)					{						Gmax2 = -G[i];						Gmax2_idx = i;					}				}				if(!is_lower_bound(i))	// d = -1				{					if(G[i] >= Gmax1)					{						Gmax1 = G[i];						Gmax1_idx = i;					}				}			}		}		if(Gmax1+Gmax2 < eps)	 		return 1;		working_set[0] = Gmax1_idx;		working_set[1] = Gmax2_idx;		return 0;	}	void do_shrinking()	{		int i,j,k;		int[] working_set = new int[2];		if(max_violating_pair(working_set)!=0) return;		i = working_set[0];		j = working_set[1];		double Gm1 = -y[j]*G[j];		double Gm2 = y[i]*G[i];		// shrink			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] >= Gm2) continue;			}			else if(is_upper_bound(k))			{				if(y[k]==+1)				{					if(G[k] >= Gm2) continue;				}				else	if(G[k] >= Gm1) 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 || -(Gm1 + Gm2) > 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] < Gm2) continue;			}			else if(is_upper_bound(k))			{				if(y[k]==+1)				{					if(G[k] < Gm2) continue;				}				else	if(G[k] < Gm1) continue;			}			else continue;			swap_index(k,active_size);			active_size++;			++k;	// look at the newcomer		}	}	double calculate_rho()	{		double r;		int nr_free = 0;		double ub = INF, lb = -INF, sum_free = 0;		for(int i=0;i<active_size;i++)		{			double yG = y[i]*G[i];			if(is_lower_bound(i))			{				if(y[i] > 0)					ub = Math.min(ub,yG);				else					lb = Math.max(lb,yG);			}			else if(is_upper_bound(i))			{				if(y[i] < 0)					ub = Math.min(ub,yG);				else					lb = Math.max(lb,yG);			}			else			{				++nr_free;				sum_free += yG;			}		}		if(nr_free>0)			r = sum_free/nr_free;		else			r = (ub+lb)/2;		return r;	}}//// Solver for nu-svm classification and regression//// additional constraint: e^T \alpha = constant//final class Solver_NU extends Solver{	private SolutionInfo si;	void Solve(int l, QMatrix Q, double[] b, byte[] y,		   double[] alpha, double Cp, double Cn, double eps,		   SolutionInfo si, int shrinking)	{		this.si = si;		super.Solve(l,Q,b,y,alpha,Cp,Cn,eps,si,shrinking);	}	// return 1 if already optimal, return 0 otherwise	int select_working_set(int[] working_set)	{		// return i,j such that y_i = y_j and		// i: maximizes -y_i * grad(f)_i, i in I_up(\alpha)		// j: minimizes the decrease of obj value		//    (if quadratic coefficeint <= 0, replace it with tau)		//    -y_j*grad(f)_j < -y_i*grad(f)_i, j in I_low(\alpha)			double Gmaxp = -INF;		double Gmaxp2 = -INF;		int Gmaxp_idx = -1;			double Gmaxn = -INF;		double Gmaxn2 = -INF;		int Gmaxn_idx = -1;			int Gmin_idx = -1;		double obj_diff_min = INF;			for(int t=0;t<active_size;t++)			if(y[t]==+1)			{				if(!is_upper_bound(t))					if(-G[t] >= Gmaxp)					{						Gmaxp = -G[t];						Gmaxp_idx = t;					}			}			else			{				if(!is_lower_bound(t))					if(G[t] >= Gmaxn)					{						Gmaxn = G[t];						Gmaxn_idx = t;					}			}			int ip = Gmaxp_idx;		int in = Gmaxn_idx;		Qfloat[] Q_ip = null;		Qfloat[] Q_in = null;		if(ip != -1) // null Q_ip not accessed: Gmaxp=-INF if ip=-1			Q_ip = Q.get_Q(ip,active_size);		if(in != -1)			Q_in = Q.get_Q(in,active_size);			for(int j=0;j<active_size;j++)		{			if(y[j]==+1)			{				if (!is_lower_bound(j))					{					double grad_diff=Gmaxp+G[j];					if (G[j] >= Gmaxp2)						Gmaxp2 = G[j];					if (grad_diff > 0)					{						double obj_diff; 						double quad_coef = Q_ip[ip]+QD[j]-2*Q_ip[j];						if (quad_coef > 0)							obj_diff = -(grad_diff*grad_diff)/quad_coef;						else							obj_diff = -(grad_diff*grad_diff)/TAU;							if (obj_diff <= obj_diff_min)						{							Gmin_idx=j;							obj_diff_min = obj_diff;						}					}				}			}			else			{				if (!is_upper_bound(j))				{					double grad_diff=Gmaxn-G[j];					if (-G[j] >= Gmaxn2)						Gmaxn2 = -G[j];					if (grad_diff > 0)					{						double obj_diff; 						double quad_coef = Q_in[in]+QD[j]-2*Q_in[j];						if (quad_coef > 0)							obj_diff = -(grad_diff*grad_diff)/quad_coef;						else							obj_diff = -(grad_diff*grad_diff)/TAU;							if (obj_diff <= obj_diff_min)						{							Gmin_idx=j;							obj_diff_min = obj_diff;						}					}				}			}		}		if(Math.max(Gmaxp+Gmaxp2,Gmaxn+Gmaxn2) < eps) 			return 1;			if(y[Gmin_idx] == +1)			working_set[0] = Gmaxp_idx;		else			working_set[0] = Gmaxn_idx;		working_set[1] = Gmin_idx;			return 0;	}	void do_shrinking()	{		double Gmax1 = -INF;	// max { -y_i * grad(f)_i | y_i = +1, i in I_up(\alpha) }		double Gmax2 = -INF;	// max { y_i * grad(f)_i | y_i = +1, i in I_low(\alpha) }		double Gmax3 = -INF;	// max { -y_i * grad(f)_i | y_i = -1, i in I_up(\alpha) }		double Gmax4 = -INF;	// max { y_i * grad(f)_i | y_i = -1, i in I_low(\alpha) } 		// find maximal violating pair first		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];			}		}		// shrinking		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 || 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		}	}		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 = Math.min(ub1,G[i]);				else if(is_upper_bound(i))					lb1 = Math.max(lb1,G[i]);				else				{					++nr_free1;					sum_free1 += G[i];				}			}			else			{				if(is_lower_bound(i))					ub2 = Math.min(ub2,G[i]);				else if(is_upper_bound(i))					lb2 = 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 extends Kernel{	private final byte[] y;	private final Cache cache;	private final Qfloat[] QD;	SVC_Q(svm_problem prob, svm_parameter param, byte[] y_)	{		super(prob.l, prob.x, param);		y = (byte[])y_.clone();		cache = new Cache(prob.l,(int)(param.cache_size*(1<<20)));		QD = new Qfloat[prob.l];		for(int i=0;i<prob.l;i++)			QD[i]= (Qfloat)kernel_function(i,i);	}	Qfloat[] get_Q(int i, int len)	{		Qfloat[][] data = new Qfloat[1][];		int start;		if((start = cache.get_data(i,data,len)) < len)		{			for(int j=start;j<len;j++)				data[0][j] = (Qfloat)(y[i]*y[j]*kernel_function(i,j));		}		return data[0];	}	Qfloat[] get_QD()	{		return QD;	}	void swap_index(int i, int j)	{		cache.swap_index(i,j);		super.swap_index(i,j);		swap(byte,y[i],y[j]);		swap(Qfloat,QD[i],QD[j]);	}}class ONE_CLASS_Q extends Kernel{	private final Cache cache;	private final Qfloat[] QD;	ONE_CLASS_Q(svm_problem prob, svm_parameter param)	{		super(prob.l, prob.x, param);		cache = new Cache(prob.l,(int)(param.cache_size*(1<<20)));		QD = new Qfloat[prob.l];		for(int i=0;i<prob.l;i++)			QD[i]= (Qfloat)kernel_function(i,i);	}	Qfloat[] get_Q(int i, int len)	{		Qfloat[][] data = new Qfloat[1][];		int start;		if((start = cache.get_data(i,data,len)) < len)		{			for(int j=start;j<len;j++)				data[0][j] = (Qfloat)kernel_function(i,j);		}		return data[0];	}	Qfloat[] get_QD()	{		return QD;	}	void swap_index(int i, int j)	{		cache.swap_index(i,j);		super.swap_index(i,j);		swap(Qfloat,QD[i],QD[j]);	}}class SVR_Q extends Kernel{	private final int l;	private final Cache cache;	private final byte[] sign;	private final int[] index;	private int next_buffer;	private Qfloat[][] buffer;	private final Qfloat[] QD;	SVR_Q(svm_problem prob, svm_parameter param)	{		super(prob.l, prob.x, param);		l = prob.l;		cache = new Cache(l,(int)(param.cache_size*(1<<20)));		QD = new Qfloat[2*l];		sign = new byte[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;			QD[k] = (Qfloat)kernel_function(k,k);			QD[k+l] = QD[k];		}		buffer = new Qfloat[2][2*l];		next_buffer = 0;	}	void swap_index(int i, int j)	{		swap(byte,sign[i],sign[j]);		swap(int,index[i],index[j]);		swap(Qfloat,QD[i],QD[j]);	}	Qfloat[] get_Q(int i, int len)	{		Qfloat[][] data = new Qfloat[1][];		int real_i = index[i];		if(cache.get_data(real_i,data,l) < l)		{			for(int j=0;j<l;j++)				data[0][j] = (Qfloat)kernel_function(real_i,j);		}		// reorder and copy		Qfloat buf[] = buffer[next_buffer];		next_buffer = 1 - next_buffer;		byte si = sign[i];		for(int j=0;j<len;j++)			buf[j] = si * sign[j] * data[0][index[j]];		return buf;	}	Qfloat[] get_QD()	{		return QD;	}}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];		byte[] y = new byte[l];		int i;		for(i=0;i<l;i++)		{			alpha[i] = 0;			minus_ones[i] = -1;			if(prob.y[i] > 0) y[i] = +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.out.print("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;		byte[] y = new byte[l];		for(i=0;i<l;i++)			if(prob.y[i]>0)				y[i] = +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] = Math.min(1.0,sum_pos);				sum_pos -= alpha[i];			}			else			{				alpha[i] = 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;