quad_coef = TAU;			double delta = (G[i]-G[j])/quad_coef;			double sum = alpha[i] + alpha[j];			alpha[i] -= delta;			alpha[j] += delta;			if(sum > C_i)			{				if(alpha[i] > C_i)				{					alpha[i] = C_i;					alpha[j] = sum - C_i;				}			}			else			{				if(alpha[j] < 0)				{					alpha[j] = 0;					alpha[i] = sum;				}			}			if(sum > C_j)			{				if(alpha[j] > C_j)				{					alpha[j] = C_j;					alpha[i] = sum - C_j;				}			}			else			{				if(alpha[i] < 0)				{					alpha[i] = 0;					alpha[j] = sum;				}			}		}		// update G		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		{			bool ui = is_upper_bound(i);			bool 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];	}	// juggle everything back	/*{		for(int i=0;i<l;i++)			while(active_set[i] != i)				swap_index(i,active_set[i]);				// or Q.swap_index(i,active_set[i]);	}*/	si->upper_bound_p = Cp;	si->upper_bound_n = Cn;	info("\noptimization finished, #iter = %d\n",iter);	delete[] b;	delete[] y;	delete[] alpha;	delete[] alpha_status;	delete[] active_set;	delete[] G;	delete[] G_bar;}// return 1 if already optimal, return 0 otherwiseint Solver::select_working_set(int &out_i, int &out_j){	// return i,j such that	// 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 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;	const Qfloat *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)/TAU;					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)/TAU;					if (obj_diff <= obj_diff_min)					{						Gmin_idx=j;						obj_diff_min = obj_diff;					}				}			}		}	}	if(Gmax+Gmax2 < eps)		return 1;	out_i = Gmax_idx;	out_j = Gmin_idx;	return 0;}// return 1 if already optimal, return 0 otherwiseint Solver::max_violating_pair(int &out_i, int &out_j){	// return i,j: maximal violating pair	double Gmax1 = -INF;		// max { -y_i * grad(f)_i | i in I_up(\alpha) }	int Gmax1_idx = -1;	double Gmax2 = -INF;		// max { y_i * grad(f)_i | i in I_low(\alpha) }	int Gmax2_idx = -1;	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;	out_i = Gmax1_idx;	out_j = Gmax2_idx;	return 0;}void Solver::do_shrinking(){	int i,j,k;	if(max_violating_pair(i,j)!=0) return;	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 Solver::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 = min(ub,yG);			else				lb = max(lb,yG);		}		else if(is_upper_bound(i))		{			if(y[i] < 0)				ub = min(ub,yG);			else				lb = 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//class Solver_NU : public Solver{public:	Solver_NU() {}	void Solve(int l, const QMatrix& Q, const double *b, const schar *y,		   double *alpha, double Cp, double Cn, double eps,		   SolutionInfo* si, int shrinking)	{		this->si = si;		Solver::Solve(l,Q,b,y,alpha,Cp,Cn,eps,si,shrinking);	}private:	SolutionInfo *si;	int select_working_set(int &i, int &j);	double calculate_rho();	void do_shrinking();};// return 1 if already optimal, return 0 otherwiseint Solver_NU::select_working_set(int &out_i, int &out_j){	// 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;	const Qfloat *Q_ip = NULL;	const 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(max(Gmaxp+Gmaxp2,Gmaxn+Gmaxn2) < eps) 		return 1;	if (y[Gmin_idx] == +1)		out_i = Gmaxp_idx;	else		out_i = Gmaxn_idx;	out_j = Gmin_idx;	return 0;}void Solver_NU::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 || 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++;