double delta_alpha_i = alpha[i] - old_alpha_i;			double delta_alpha_j = alpha[j] - old_alpha_j;			for(int k=0;k<active_size;k++)			{				G[k] += Q_i[k]*delta_alpha_i + Q_j[k]*delta_alpha_j;			}			// update alpha_status and G_bar			{				boolean ui = is_upper_bound(i);				boolean uj = is_upper_bound(j);				update_alpha_status(i);				update_alpha_status(j);				int k;				if(ui != is_upper_bound(i))				{					Q_i = Q.get_Q(i,l);					if(ui)						for(k=0;k<l;k++)							G_bar[k] -= C_i * Q_i[k];					else						for(k=0;k<l;k++)							G_bar[k] += C_i * Q_i[k];				}				if(uj != is_upper_bound(j))				{					Q_j = Q.get_Q(j,l);					if(uj)						for(k=0;k<l;k++)							G_bar[k] -= C_j * Q_j[k];					else						for(k=0;k<l;k++)							G_bar[k] += C_j * Q_j[k];				}			}		}		// calculate rho		si.rho = calculate_rho();		// calculate objective value		{			double v = 0;			int i;			for(i=0;i<l;i++)				v += alpha[i] * (G[i] + b[i]);			si.obj = v/2;		}		// put back the solution		{			for(int i=0;i<l;i++)				alpha_[active_set[i]] = alpha[i];		}		si.upper_bound_p = Cp;		si.upper_bound_n = Cn;		System.out.print("\noptimization finished, #iter = "+iter+"\n");	}	// return 1 if already optimal, return 0 otherwise	int select_working_set(int[] working_set)	{		// return i,j such that		// i: maximizes -y_i * grad(f)_i, i in I_up(\alpha)		// j: mimimizes 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 Gmax = -INF;		double Gmax2 = -INF;		int Gmax_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] >= Gmax)					{						Gmax = -G[t];						Gmax_idx = t;					}			}			else			{				if(!is_lower_bound(t))					if(G[t] >= Gmax)					{						Gmax = G[t];						Gmax_idx = t;					}			}			int i = Gmax_idx;		float[] Q_i = null;		if(i != -1) // null Q_i not accessed: Gmax=-INF if i=-1			Q_i = Q.get_Q(i,active_size);			for(int j=0;j<active_size;j++)		{			if(y[j]==+1)			{				if (!is_lower_bound(j))				{					double grad_diff=Gmax+G[j];					if (G[j] >= Gmax2)						Gmax2 = G[j];					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)/1e-12;							if (obj_diff <= obj_diff_min)						{							Gmin_idx=j;							obj_diff_min = obj_diff;						}					}				}			}			else			{				if (!is_upper_bound(j))				{					double grad_diff= Gmax-G[j];					if (-G[j] >= Gmax2)						Gmax2 = -G[j];					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)/1e-12;							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;		float[] Q_ip = null;		float[] 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)/1e-12;							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)/1e-12;							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;