++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)));		QD = new Qfloat[prob.l];		for(int i=0;i<prob.l;i++)			QD[i]= (Qfloat)(this->*kernel_function)(i,i);	}		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;	}	Qfloat *get_QD() const	{		return QD;	}	void swap_index(int i, int j) const	{		cache->swap_index(i,j);		Kernel::swap_index(i,j);		swap(y[i],y[j]);		swap(QD[i],QD[j]);	}	~SVC_Q()	{		delete[] y;		delete cache;		delete[] QD;	}private:	schar *y;	Cache *cache;	Qfloat *QD;};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)));		QD = new Qfloat[prob.l];		for(int i=0;i<prob.l;i++)			QD[i]= (Qfloat)(this->*kernel_function)(i,i);	}		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;	}	Qfloat *get_QD() const	{		return QD;	}	void swap_index(int i, int j) const	{		cache->swap_index(i,j);		Kernel::swap_index(i,j);		swap(QD[i],QD[j]);	}	~ONE_CLASS_Q()	{		delete cache;		delete[] QD;	}private:	Cache *cache;	Qfloat *QD;};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)));		QD = new Qfloat[2*l];		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;			QD[k]= (Qfloat)(this->*kernel_function)(k,k);			QD[k+l]=QD[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]);		swap(QD[i],QD[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;	}	Qfloat *get_QD() const	{		return QD;	}	~SVR_Q()	{		delete cache;		delete[] sign;		delete[] index;		delete[] buffer[0];		delete[] buffer[1];		delete[] QD;	}private:	int l;	Cache *cache;	schar *sign;	int *index;	mutable int next_buffer;	Qfloat *buffer[2];	Qfloat *QD;};//// 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;	if(n<prob->l)		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;		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;	s.Solve(2*l, 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 += fabs(alpha[i]);	}	info("nu = %f\n",sum_alpha/(param->C*l));	delete[] alpha2;	delete[] linear_term;	delete[] y;}static void solve_nu_svr(	const svm_problem *prob, const 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];	schar *y = new schar[2*l];	int i;	double sum = C * param->nu * l / 2;	for(i=0;i<l;i++)	{		alpha2[i] = alpha2[i+l] = min(sum,C);		sum -= alpha2[i];		linear_term[i] = - prob->y[i];		y[i] = 1;		linear_term[i+l] = prob->y[i];		y[i+l] = -1;	}	Solver_NU s;	s.Solve(2*l, SVR_Q(*prob,*param), linear_term, y,		alpha2, C, C, param->eps, si, param->shrinking);	info("epsilon = %f\n",-si->r);	for(i=0;i<l;i++)		alpha[i] = alpha2[i] - alpha2[i+l];	delete[] alpha2;	delete[] linear_term;	delete[] y;}//// decision_function//struct decision_function{	double *alpha;	double rho;	};decision_function svm_train_one(	const svm_problem *prob, const svm_parameter *param,	double Cp, double Cn){	double *alpha = Malloc(double,prob->l);	Solver::SolutionInfo si;	switch(param->svm_type)	{		case C_SVC:			solve_c_svc(prob,param,alpha,&si,Cp,Cn);			break;		case NU_SVC:			solve_nu_svc(prob,param,alpha,&si);			break;		case ONE_CLASS:			solve_one_class(prob,param,alpha,&si);			break;		case EPSILON_SVR:			solve_epsilon_svr(prob,param,alpha,&si);			break;		case NU_SVR:			solve_nu_svr(prob,param,alpha,&si);			break;	}	info("obj = %f, rho = %f\n",si.obj,si.rho);	// output SVs	int nSV = 0;	int nBSV = 0;	for(int i=0;i<prob->l;i++)	{		if(fabs(alpha[i]) > 0)		{			++nSV;			if(prob->y[i] > 0)			{				if(fabs(alpha[i]) >= si.upper_bound_p)					++nBSV;			}			else			{				if(fabs(alpha[i]) >= si.upper_bound_n)					++nBSV;			}		}	}	info("nSV = %d, nBSV = %d\n",nSV,nBSV);	decision_function f;	f.alpha = alpha;	f.rho = si.rho;	return f;}//// svm_model//struct svm_model{	svm_parameter param;	// parameter	int nr_class;		// number of classes, = 2 in regression/one class svm	int l;			// total #SV	svm_node **SV;		// SVs (SV[l])	double **sv_coef;	// coefficients for SVs in decision functions (sv_coef[n-1][l])	double *rho;		// constants in decision functions (rho[n*(n-1)/2])	double *probA;          // pariwise probability information	double *probB;	// for classification only	int *label;		// label of each class (label[n])	int *nSV;		// number of SVs for each class (nSV[n])				// nSV[0] + nSV[1] + ... + nSV[n-1] = l	// XXX	int free_sv;		// 1 if svm_model is created by svm_load_model				// 0 if svm_model is created by svm_train};// Platt's binary SVM Probablistic Output: an improvement from Lin et al.void sigmoid_train(	int l, const double *dec_values, const double *labels, 	double& A, double& B){	double prior1=0, prior0 = 0;	int i;	for (i=0;i<l;i++)		if (labels[i] > 0) prior1+=1;		else prior0+=1;		int max_iter=100; 	// Maximal number of iterations	double min_step=1e-10;	// Minimal step taken in line search	double sigma=1e-3;	// For numerically strict PD of Hessian	double eps=1e-5;	double hiTarget=(prior1+1.0)/(prior1+2.0);	double loTarget=1/(prior0+2.0);	double *t=Malloc(double,l);	double fApB,p,q,h11,h22,h21,g1,g2,det,dA,dB,gd,stepsize;	double newA,newB,newf,d1,d2;	int iter; 		// Initial Point and Initial Fun Value	A=0.0; B=log((prior0+1.0)/(prior1+1.0));	double fval = 0.0;	for (i=0;i<l;i++)	{		if (labels[i]>0) t[i]=hiTarget;		else t[i]=loTarget;		fApB = dec_values[i]*A+B;		if (fApB>=0)			fval += t[i]*fApB + log(1+exp(-fApB));		else			fval += (t[i] - 1)*fApB +log(1+exp(fApB));	}	for (iter=0;iter<max_iter;iter++)	{		// Update Gradient and Hessian (use H' = H + sigma I)		h11=sigma; // numerically ensures strict PD		h22=sigma;		h21=0.0;g1=0.0;g2=0.0;		for (i=0;i<l;i++)		{			fApB = dec_values[i]*A+B;			if (fApB >= 0)			{				p=exp(-fApB)/(1.0+exp(-fApB));				q=1.0/(1.0+exp(-fApB));			}			else			{				p=1.0/(1.0+exp(fApB));				q=exp(fApB)/(1.0+exp(fApB));			}			d2=p*q;			h11+=dec_values[i]*dec_values[i]*d2;			h22+=d2;			h21+=dec_values[i]*d2;			d1=t[i]-p;			g1+=dec_values[i]*d1;			g2+=d1;		}		// Stopping Criteria		if (fabs(g1)<eps && fabs(g2)<eps)			break;		// Finding Newton direction: -inv(H') * g		det=h11*h22-h21*h21;		dA=-(h22*g1 - h21 * g2) / det;		dB=-(-h21*g1+ h11 * g2) / det;		gd=g1*dA+g2*dB;		stepsize = 1; 		// Line Search		while (stepsize >= min_step)		{			newA = A + stepsize * dA;			newB = B + stepsize * dB;			// New function value			newf = 0.0;			for (i=0;i<l;i++)			{				fApB = dec_values[i]*newA+newB;				if (fApB >= 0)					newf += t[i]*fApB + log(1+exp(-fApB));				else					newf += (t[i] - 1)*fApB +log(1+exp(fApB));			}			// Check sufficient decrease			if (newf<fval+0.0001*stepsize*gd)			{				A=newA;B=newB;fval=newf;				break;