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 which maximize -grad(f)^T d , under constraint	// if alpha_i == C, d != +1	// if alpha_i == 0, d != -1	double Gmax1 = -INF;		// max { -grad(f)_i * d | y_i*d = +1 }	int Gmax1_idx = -1;	double Gmax2 = -INF;		// max { -grad(f)_i * d | y_i*d = -1 }	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(select_working_set(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 Kernel& 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();};int Solver_NU::select_working_set(int &out_i, int &out_j){	// 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 { -grad(f)_i * d | y_i = +1, d = +1 }	int Gmax1_idx = -1;	double Gmax2 = -INF;	// max { -grad(f)_i * d | y_i = +1, d = -1 }	int Gmax2_idx = -1;	double Gmax3 = -INF;	// max { -grad(f)_i * d | y_i = -1, d = +1 }	int Gmax3_idx = -1;	double Gmax4 = -INF;	// max { -grad(f)_i * d | y_i = -1, d = -1 }	int Gmax4_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] > 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(max(Gmax1+Gmax2,Gmax3+Gmax4) < eps) 		return 1;	if(Gmax1+Gmax2 > Gmax3+Gmax4)	{		out_i = Gmax1_idx;		out_j = Gmax2_idx;	}	else	{		out_i = Gmax3_idx;		out_j = Gmax4_idx;	}	return 0;}void Solver_NU::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 || 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 Solver_NU::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 = min(ub1,G[i]);			else if(is_upper_bound(i))				lb1 = max(lb1,G[i]);			else			{				++nr_free1;				sum_free1 += G[i];			}		}		else		{			if(is_lower_bound(i))				ub2 = min(ub2,G[i]);			else if(is_upper_bound(i))				lb2 = 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: public Kernel{ public:	SVC_Q(const svm_problem& prob, const svm_parameter& param, const schar *y_)	:Kernel(prob.l, prob.x, param)	{		clone(y,y_,prob.l);		cache = new Cache(prob.l,(int)(param.cache_size*(1<<20)));	}		Qfloat *get_Q(int i, int len) const	{		Qfloat *data;		int start;		if((start = cache->get_data(i,&data,len)) < len)		{			for(int j=start;j<len;j++)				data[j] = (Qfloat)(y[i]*y[j]*(this->*kernel_function)(i,j));		}		return data;	}	void swap_index(int i, int j) const	{		cache->swap_index(i,j);		Kernel::swap_index(i,j);		swap(y[i],y[j]);	}	~SVC_Q()	{		delete[] y;		delete cache;	}private:	schar *y;	Cache *cache;};class ONE_CLASS_Q: public Kernel{public:	ONE_CLASS_Q(const svm_problem& prob, const svm_parameter& param)	:Kernel(prob.l, prob.x, param)	{		cache = new Cache(prob.l,(int)(param.cache_size*(1<<20)));	}		Qfloat *get_Q(int i, int len) const	{		Qfloat *data;		int start;		if((start = cache->get_data(i,&data,len)) < len)		{			for(int j=start;j<len;j++)				data[j] = (Qfloat)(this->*kernel_function)(i,j);		}		return data;	}	void swap_index(int i, int j) const	{		cache->swap_index(i,j);		Kernel::swap_index(i,j);	}	~ONE_CLASS_Q()	{		delete cache;	}private:	Cache *cache;};class SVR_Q: public Kernel{ public:	SVR_Q(const svm_problem& prob, const svm_parameter& param)	:Kernel(prob.l, prob.x, param)	{		l = prob.l;		cache = new Cache(l,(int)(param.cache_size*(1<<20)));		sign = new schar[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[0] = new Qfloat[2*l];		buffer[1] = new Qfloat[2*l];		next_buffer = 0;	}	void swap_index(int i, int j) const	{		swap(sign[i],sign[j]);		swap(index[i],index[j]);	}		Qfloat *get_Q(int i, int len) const	{		Qfloat *data;		int real_i = index[i];		if(cache->get_data(real_i,&data,l) < l)		{			for(int j=0;j<l;j++)				data[j] = (Qfloat)(this->*kernel_function)(real_i,j);		}		// reorder and copy		Qfloat *buf = buffer[next_buffer];		next_buffer = 1 - next_buffer;		schar si = sign[i];		for(int j=0;j<len;j++)			buf[j] = si * sign[j] * data[index[j]];		return buf;	}	~SVR_Q()	{		delete cache;		delete[] sign;		delete[] index;		delete[] buffer[0];		delete[] buffer[1];	}private:	int l;	Cache *cache;	schar *sign;	int *index;	mutable int next_buffer;	Qfloat* buffer[2];};//// construct and solve various formulations//static void solve_c_svc(	const svm_problem *prob, const svm_parameter* param,	double *alpha, Solver::SolutionInfo* si, double Cp, double Cn){	int l = prob->l;	double *minus_ones = new double[l];	schar *y = new schar[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;	s.Solve(l, 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)		info("nu = %f\n", sum_alpha/(Cp*prob->l));	for(i=0;i<l;i++)		alpha[i] *= y[i];	delete[] minus_ones;	delete[] y;}static void solve_nu_svc(	const svm_problem *prob, const svm_parameter *param,	double *alpha, Solver::SolutionInfo* si){	int i;	int l = prob->l;	double nu = param->nu;	schar *y = new schar[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] = min(1.0,sum_pos);			sum_pos -= alpha[i];		}		else		{			alpha[i] = 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;	s.Solve(l, SVC_Q(*prob,*param,y), zeros, y,		alpha, 1.0, 1.0, param->eps, si,  param->shrinking);	double r = si->r;	info("C = %f\n",1/r);	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;	delete[] y;	delete[] zeros;}static void solve_one_class(	const svm_problem *prob, const svm_parameter *param,	double *alpha, Solver::SolutionInfo* si){	int l = prob->l;	double *zeros = new double[l];	schar *ones = new schar[l];	int i;	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;	s.Solve(l, ONE_CLASS_Q(*prob,*param), zeros, ones,		alpha, 1.0, 1.0, param->eps, si, param->shrinking);	delete[] zeros;	delete[] ones;}static void solve_epsilon_svr(	const svm_problem *prob, const 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];	schar *y = new schar[2*l];	int i;	for(i=0;i<l;i++)	{		alpha2[i] = 0;