text: #ifdef GPU
  #include <cuda_runtime.h>
  #include <cublas_v2.h>
  #include <helper_cuda.h>
  #include <iostream>
  #include <assert.h>
  extern cublasHandle_t handle;
  #endif // GPU
  
  #include <math.h>
  #include <stdio.h>
  #include <stdlib.h>
  #include <string.h>
  #include <stdarg.h>
  #include <locale.h>
  #include "linear.h"
  #include "tron.h"
  typedef signed char schar;
  template <class T> static inline void swap(T& x, T& y) { T t=x; x=y; y=t; }
  #ifndef min
  template <class T> static inline T min(T x,T y) { return (x<y)?x:y; }
  #endif
  #ifndef max
  template <class T> static inline T max(T x,T y) { return (x>y)?x:y; }
  #endif
  template <class S, class T> static inline void clone(T*& dst, S* src, int n)
  {
  	dst = new T[n];
  	memcpy((void *)dst,(void *)src,sizeof(T)*n);
  }
  #define Malloc(type,n) (type *)malloc((n)*sizeof(type))
  #define INF HUGE_VAL
  
  static void print_string_stdout(const char *s)
  {
  	fputs(s,stdout);
  	fflush(stdout);
  }
  
  static void (*liblinear_print_string) (const char *) = &print_string_stdout;
  
  #if 1
  static void info(const char *fmt,...)
  {
  	char buf[BUFSIZ];
  	va_list ap;
  	va_start(ap,fmt);
  	vsprintf(buf,fmt,ap);
  	va_end(ap);
  	(*liblinear_print_string)(buf);
  }
  #else
  static void info(const char *fmt,...) {}
  #endif
  
  class l2r_lr_fun: public function
  {
  public:
  	l2r_lr_fun(const problem *prob, double *C);
  	~l2r_lr_fun();
  
  	double fun(double *w);
  	void grad(double *w, double *g);
  	void Hv(double *s, double *Hs);
  
  	int get_nr_variable(void);
  
  private:
  	void Xv(double *v, double *Xv);
  	void XTv(double *v, double *XTv);
  
  	double *C;
  	double *z;
  	double *D;
  	const problem *prob;
  
  #ifdef GPU
  	static double *X_dev;
  	static double *v_dev;
  	static double *Xv_dev;
  #endif // GPU
  };
  
  #ifdef GPU
  double* l2r_lr_fun::X_dev = NULL;
  double* l2r_lr_fun::v_dev = NULL;
  double* l2r_lr_fun::Xv_dev = NULL;
  #endif // GPU
  
  l2r_lr_fun::l2r_lr_fun(const problem *prob, double *C)
  {
  	int l=prob->l;
  
  	this->prob = prob;
  
  	z = new double[l];
  	D = new double[l];
  	this->C = C;
  
  #ifdef GPU
  	if(X_dev == NULL) {
  		printf("Copying (and allocating) the data matrix to GPU...\n");
  		double gb = 1.0*sizeof(double)*prob->n*prob->l/1024/1024/1024;
  		printf("Size of data matrix is %d bytes (i.e. %lf GB)\n", sizeof(double)*prob->n*prob->l, gb);
  
  		double* X = new double [prob->l * prob->n];
  
  		for(int i = 0; i < prob->l; ++i)
  			for(int j = 0; (prob->x)[i][j].index >=0; ++j)
  				X[i + ((prob->x)[i][j].index-1) * prob->l] = (prob->x)[i][j].value; 
  
  		assert(cudaMalloc((void **)&X_dev, prob->n*prob->l * sizeof(X_dev[0])) == cudaSuccess);
  		assert(cudaMalloc((void **)&v_dev, max(prob->n, prob->l) * sizeof(v_dev[0])) == cudaSuccess);
  		assert(cudaMalloc((void **)&Xv_dev, max(prob->n, prob->l) * sizeof(Xv_dev[0])) == cudaSuccess);
  
  	    assert(cublasSetVector(prob->n*prob->l, sizeof(X_dev[0]), X, 1, X_dev, 1) == CUBLAS_STATUS_SUCCESS);
  		printf("GPU allocation and data transfer finished.\n");
  	}
  #endif // GPU
  }
  
  l2r_lr_fun::~l2r_lr_fun()
  {
  	delete[] z;
  	delete[] D;
  }
  
  
  double l2r_lr_fun::fun(double *w)
  {
  	int i;
  	double f=0;
  	double *y=prob->y;
  	int l=prob->l;
  	int w_size=get_nr_variable();
  
  	Xv(w, z);
  
  	for(i=0;i<w_size;i++)
  		f += w[i]*w[i];
  	f /= 2.0;
  	for(i=0;i<l;i++)
  	{
  		double yz = y[i]*z[i];
  		if (yz >= 0)
  			f += C[i]*log(1 + exp(-yz));
  		else
  			f += C[i]*(-yz+log(1 + exp(yz)));
  	}
  
  	return(f);
  }
  
  void l2r_lr_fun::grad(double *w, double *g)
  {
  	int i;
  	double *y=prob->y;
  	int l=prob->l;
  	int w_size=get_nr_variable();
  
  	for(i=0;i<l;i++)
  	{
  		z[i] = 1/(1 + exp(-y[i]*z[i]));
  		D[i] = z[i]*(1-z[i]);
  		z[i] = C[i]*(z[i]-1)*y[i];
  	}
  	XTv(z, g);
  
  	for(i=0;i<w_size;i++)
  		g[i] = w[i] + g[i];
  }
  
  int l2r_lr_fun::get_nr_variable(void)
  {
  	return prob->n;
  }
  
  void l2r_lr_fun::Hv(double *s, double *Hs)
  {
  	int i;
  	int l=prob->l;
  	int w_size=get_nr_variable();
  	double *wa = new double[l];
  
  	Xv(s, wa);
  	for(i=0;i<l;i++)
  		wa[i] = C[i]*D[i]*wa[i];
  
  	XTv(wa, Hs);
  	for(i=0;i<w_size;i++)
  		Hs[i] = s[i] + Hs[i];
  	delete[] wa;
  }
  
  #ifdef GPU
  void l2r_lr_fun::Xv(double *v, double *Xv)
  {
  	int l = prob->l;
  	int n = prob->n;
  	int lda = l;
  	int inc = 1;
  	double beta = 0.0f;
  	double alpha = 1.0f;
  
      assert(cublasSetVector(n, sizeof(v_dev[0]), v, 1, v_dev, 1) == CUBLAS_STATUS_SUCCESS);
  	assert(cublasDgemv(handle, CUBLAS_OP_N, l, n, &alpha, X_dev, lda, v_dev, inc, &beta, Xv_dev, inc) == CUBLAS_STATUS_SUCCESS);
      assert(cublasGetVector(l, sizeof(Xv[0]), Xv_dev, 1, Xv, 1) == CUBLAS_STATUS_SUCCESS);
  }
  
  #else
  
  void l2r_lr_fun::Xv(double *v, double *Xv)
  {
  	int i;
  	int l=prob->l;
  	feature_node **x=prob->x;
  
  	for(i=0;i<l;i++)
  	{
  		feature_node *s=x[i];
  		Xv[i]=0;
  		while(s->index!=-1)
  		{
  			Xv[i]+=v[s->index-1]*s->value;
  			s++;
  		}
  	}
  }
  #endif // GPU
  
  #ifdef GPU
  
  void l2r_lr_fun::XTv(double *v, double *XTv)
  {
  	int l=prob->l;
  	int n=prob->n;
  
  	int lda = l;
  	int inc = 1;
  	double beta = 0.0f;
  	double alpha = 1.0f;
  
      assert(cublasSetVector(l, sizeof(v_dev[0]), v, 1, v_dev, 1) == CUBLAS_STATUS_SUCCESS);
  	assert(cublasDgemv(handle, CUBLAS_OP_T, l, n, &alpha, X_dev, lda, v_dev, inc, &beta, Xv_dev, inc) == CUBLAS_STATUS_SUCCESS);
      assert(cublasGetVector(n, sizeof(XTv[0]), Xv_dev, 1, XTv, 1) == CUBLAS_STATUS_SUCCESS);
  }
  
  #else
  
  void l2r_lr_fun::XTv(double *v, double *XTv)
  {
  	int i;
  	int l=prob->l;
  	int w_size=get_nr_variable();
  	feature_node **x=prob->x;
  
  	for(i=0;i<w_size;i++)
  		XTv[i]=0;
  	for(i=0;i<l;i++)
  	{
  		feature_node *s=x[i];
  		while(s->index!=-1)
  		{
  			XTv[s->index-1]+=v[i]*s->value;
  			s++;
  		}
  	}
  }
  #endif // GPU
  
  class l2r_l2_svc_fun: public function
  {
  public:
  	l2r_l2_svc_fun(const problem *prob, double *C);
  	~l2r_l2_svc_fun();
  
  	double fun(double *w);
  	void grad(double *w, double *g);
  	void Hv(double *s, double *Hs);
  
  	int get_nr_variable(void);
  
  protected:
  	void Xv(double *v, double *Xv);
  	void subXv(double *v, double *Xv);
  	void subXTv(double *v, double *XTv);
  
  	double *C;
  	double *z;
  	double *D;
  	int *I;
  	int sizeI;
  	const problem *prob;
  };
  
  l2r_l2_svc_fun::l2r_l2_svc_fun(const problem *prob, double *C)
  {
  	int l=prob->l;
  
  	this->prob = prob;
  
  	z = new double[l];
  	D = new double[l];
  	I = new int[l];
  	this->C = C;
  }
  
  l2r_l2_svc_fun::~l2r_l2_svc_fun()
  {
  	delete[] z;
  	delete[] D;
  	delete[] I;
  }
  
  double l2r_l2_svc_fun::fun(double *w)
  {
  	int i;
  	double f=0;
  	double *y=prob->y;
  	int l=prob->l;
  	int w_size=get_nr_variable();
  
  	Xv(w, z);
  
  	for(i=0;i<w_size;i++)
  		f += w[i]*w[i];
  	f /= 2.0;
  	for(i=0;i<l;i++)
  	{
  		z[i] = y[i]*z[i];
  		double d = 1-z[i];
  		if (d > 0)
  			f += C[i]*d*d;
  	}
  
  	return(f);
  }
  
  void l2r_l2_svc_fun::grad(double *w, double *g)
  {
  	int i;
  	double *y=prob->y;
  	int l=prob->l;
  	int w_size=get_nr_variable();
  
  	sizeI = 0;
  	for (i=0;i<l;i++)
  		if (z[i] < 1)
  		{
  			z[sizeI] = C[i]*y[i]*(z[i]-1);
  			I[sizeI] = i;
  			sizeI++;
  		}
  	subXTv(z, g);
  
  	for(i=0;i<w_size;i++)
  		g[i] = w[i] + 2*g[i];
  }
  
  int l2r_l2_svc_fun::get_nr_variable(void)
  {
  	return prob->n;
  }
  
  void l2r_l2_svc_fun::Hv(double *s, double *Hs)
  {
  	int i;
  	int w_size=get_nr_variable();
  	double *wa = new double[sizeI];
  
  	subXv(s, wa);
  	for(i=0;i<sizeI;i++)
  		wa[i] = C[I[i]]*wa[i];
  
  	subXTv(wa, Hs);
  	for(i=0;i<w_size;i++)
  		Hs[i] = s[i] + 2*Hs[i];
  	delete[] wa;
  }
  
  void l2r_l2_svc_fun::Xv(double *v, double *Xv)
  {
  	int i;
  	int l=prob->l;
  	feature_node **x=prob->x;
  
  	for(i=0;i<l;i++)
  	{
  		feature_node *s=x[i];
  		Xv[i]=0;
  		while(s->index!=-1)
  		{
  			Xv[i]+=v[s->index-1]*s->value;
  			s++;
  		}
  	}
  }
  
  void l2r_l2_svc_fun::subXv(double *v, double *Xv)
  {
  	int i;
  	feature_node **x=prob->x;
  
  	for(i=0;i<sizeI;i++)
  	{
  		feature_node *s=x[I[i]];
  		Xv[i]=0;
  		while(s->index!=-1)
  		{
  			Xv[i]+=v[s->index-1]*s->value;
  			s++;
  		}
  	}
  }
  
  void l2r_l2_svc_fun::subXTv(double *v, double *XTv)
  {
  	int i;
  	int w_size=get_nr_variable();
  	feature_node **x=prob->x;
  
  	for(i=0;i<w_size;i++)
  		XTv[i]=0;
  	for(i=0;i<sizeI;i++)
  	{
  		feature_node *s=x[I[i]];
  		while(s->index!=-1)
  		{
  			XTv[s->index-1]+=v[i]*s->value;
  			s++;
  		}
  	}
  }
  
  class l2r_l2_svr_fun: public l2r_l2_svc_fun
  {
  public:
  	l2r_l2_svr_fun(const problem *prob, double *C, double p);
  
  	double fun(double *w);
  	void grad(double *w, double *g);
  
  private:
  	double p;
  };
  
  l2r_l2_svr_fun::l2r_l2_svr_fun(const problem *prob, double *C, double p):
  	l2r_l2_svc_fun(prob, C)
  {
  	this->p = p;
  }
  
  double l2r_l2_svr_fun::fun(double *w)
  {
  	int i;
  	double f=0;
  	double *y=prob->y;
  	int l=prob->l;
  	int w_size=get_nr_variable();
  	double d;
  
  	Xv(w, z);
  
  	for(i=0;i<w_size;i++)
  		f += w[i]*w[i];
  	f /= 2;
  	for(i=0;i<l;i++)
  	{
  		d = z[i] - y[i];
  		if(d < -p)
  			f += C[i]*(d+p)*(d+p);
  		else if(d > p)
  			f += C[i]*(d-p)*(d-p);
  	}
  
  	return(f);
  }
  
  void l2r_l2_svr_fun::grad(double *w, double *g)
  {
  	int i;
  	double *y=prob->y;
  	int l=prob->l;
  	int w_size=get_nr_variable();
  	double d;
  
  	sizeI = 0;
  	for(i=0;i<l;i++)
  	{
  		d = z[i] - y[i];
  
  		// generate index set I
  		if(d < -p)
  		{
  			z[sizeI] = C[i]*(d+p);
  			I[sizeI] = i;
  			sizeI++;
  		}
  		else if(d > p)
  		{
  			z[sizeI] = C[i]*(d-p);
  			I[sizeI] = i;
  			sizeI++;
  		}
  
  	}
  	subXTv(z, g);
  
  	for(i=0;i<w_size;i++)
  		g[i] = w[i] + 2*g[i];
  }
  
  // A coordinate descent algorithm for 
  // multi-class support vector machines by Crammer and Singer
  //
  //  min_{\alpha}  0.5 \sum_m ||w_m(\alpha)||^2 + \sum_i \sum_m e^m_i alpha^m_i
  //    s.t.     \alpha^m_i <= C^m_i \forall m,i , \sum_m \alpha^m_i=0 \forall i
  // 
  //  where e^m_i = 0 if y_i  = m,
  //        e^m_i = 1 if y_i != m,
  //  C^m_i = C if m  = y_i, 
  //  C^m_i = 0 if m != y_i, 
  //  and w_m(\alpha) = \sum_i \alpha^m_i x_i 
  //
  // Given: 
  // x, y, C
  // eps is the stopping tolerance
  //
  // solution will be put in w
  //
  // See Appendix of LIBLINEAR paper, Fan et al. (2008)
  
  #define GETI(i) ((int) prob->y[i])
  // To support weights for instances, use GETI(i) (i)
  
  class Solver_MCSVM_CS
  {
  	public:
  		Solver_MCSVM_CS(const problem *prob, int nr_class, double *C, double eps=0.1, int max_iter=100000);
  		~Solver_MCSVM_CS();
  		void Solve(double *w);
  	private:
  		void solve_sub_problem(double A_i, int yi, double C_yi, int active_i, double *alpha_new);
  		bool be_shrunk(int i, int m, int yi, double alpha_i, double minG);
  		double *B, *C, *G;
  		int w_size, l;
  		int nr_class;
  		int max_iter;
  		double eps;
  		const problem *prob;
  };
  
  Solver_MCSVM_CS::Solver_MCSVM_CS(const problem *prob, int nr_class, double *weighted_C, double eps, int max_iter)
  {
  	this->w_size = prob->n;
  	this->l = prob->l;
  	this->nr_class = nr_class;
  	this->eps = eps;
  	this->max_iter = max_iter;
  	this->prob = prob;
  	this->B = new double[nr_class];
  	this->G = new double[nr_class];
  	this->C = weighted_C;
  }
  
  Solver_MCSVM_CS::~Solver_MCSVM_CS()
  {
  	delete[] B;
  	delete[] G;
  }
  
  int compare_double(const void *a, const void *b)
  {
  	if(*(double *)a > *(double *)b)
  		return -1;
  	if(*(double *)a < *(double *)b)
  		return 1;
  	return 0;
  }
  
  void Solver_MCSVM_CS::solve_sub_problem(double A_i, int yi, double C_yi, int active_i, double *alpha_new)
  {
  	int r;
  	double *D;
  
  	clone(D, B, active_i);
  	if(yi < active_i)
  		D[yi] += A_i*C_yi;
  	qsort(D, active_i, sizeof(double), compare_double);
  
  	double beta = D[0] - A_i*C_yi;
  	for(r=1;r<active_i && beta<r*D[r];r++)
  		beta += D[r];
  	beta /= r;
  
  	for(r=0;r<active_i;r++)
  	{
  		if(r == yi)
  			alpha_new[r] = min(C_yi, (beta-B[r])/A_i);
  		else
  			alpha_new[r] = min((double)0, (beta - B[r])/A_i);
  	}
  	delete[] D;
  }
  
  bool Solver_MCSVM_CS::be_shrunk(int i, int m, int yi, double alpha_i, double minG)
  {
  	double bound = 0;
  	if(m == yi)
  		bound = C[GETI(i)];
  	if(alpha_i == bound && G[m] < minG)
  		return true;
  	return false;
  }
  
  void Solver_MCSVM_CS::Solve(double *w)
  {
  	int i, m, s;
  	int iter = 0;
  	double *alpha =  new double[l*nr_class];
  	double *alpha_new = new double[nr_class];
  	int *index = new int[l];
  	double *QD = new double[l];
  	int *d_ind = new int[nr_class];
  	double *d_val = new double[nr_class];
  	int *alpha_index = new int[nr_class*l];
  	int *y_index = new int[l];
  	int active_size = l;
  	int *active_size_i = new int[l];
  	double eps_shrink = max(10.0*eps, 1.0); // stopping tolerance for shrinking
  	bool start_from_all = true;
  
  	// Initial alpha can be set here. Note that 
  	// sum_m alpha[i*nr_class+m] = 0, for all i=1,...,l-1
  	// alpha[i*nr_class+m] <= C[GETI(i)] if prob->y[i] == m
  	// alpha[i*nr_class+m] <= 0 if prob->y[i] != m
  	// If initial alpha isn't zero, uncomment the for loop below to initialize w
  	for(i=0;i<l*nr_class;i++)
  		alpha[i] = 0;
  
  	for(i=0;i<w_size*nr_class;i++)
  		w[i] = 0;
  	for(i=0;i<l;i++)
  	{
  		for(m=0;m<nr_class;m++)
  			alpha_index[i*nr_class+m] = m;
  		feature_node *xi = prob->x[i];
  		QD[i] = 0;
  		while(xi->index != -1)
  		{
  			double val = xi->value;
  			QD[i] += val*val;
  
  			// Uncomment the for loop if initial alpha isn't zero
  			// for(m=0; m<nr_class; m++)
  			//	w[(xi->index-1)*nr_class+m] += alpha[i*nr_class+m]*val;
  			xi++;
  		}
  		active_size_i[i] = nr_class;
  		y_index[i] = (int)prob->y[i];
  		index[i] = i;
  	}
  
  	while(iter < max_iter)
  	{
  		double stopping = -INF;
  		for(i=0;i<active_size;i++)
  		{
  			int j = i+rand()%(active_size-i);
  			swap(index[i], index[j]);
  		}
  		for(s=0;s<active_size;s++)
  		{
  			i = index[s];
  			double Ai = QD[i];
  			double *alpha_i = &alpha[i*nr_class];
  			int *alpha_index_i = &alpha_index[i*nr_class];
  
  			if(Ai > 0)
  			{
  				for(m=0;m<active_size_i[i];m++)
  					G[m] = 1;
  				if(y_index[i] < active_size_i[i])
  					G[y_index[i]] = 0;
  
  				feature_node *xi = prob->x[i];
  				while(xi->index!= -1)
  				{
  					double *w_i = &w[(xi->index-1)*nr_class];
  					for(m=0;m<active_size_i[i];m++)
  						G[m] += w_i[alpha_index_i[m]]*(xi->value);
  					xi++;
  				}
  
  				double minG = INF;
  				double maxG = -INF;
  				for(m=0;m<active_size_i[i];m++)
  				{
  					if(alpha_i[alpha_index_i[m]] < 0 && G[m] < minG)
  						minG = G[m];
  					if(G[m] > maxG)
  						maxG = G[m];
  				}
  				if(y_index[i] < active_size_i[i])
  					if(alpha_i[(int) prob->y[i]] < C[GETI(i)] && G[y_index[i]] < minG)
  						minG = G[y_index[i]];
  
  				for(m=0;m<active_size_i[i];m++)
  				{
  					if(be_shrunk(i, m, y_index[i], alpha_i[alpha_index_i[m]], minG))
  					{
  						active_size_i[i]--;
  						while(active_size_i[i]>m)
  						{
  							if(!be_shrunk(i, active_size_i[i], y_index[i],
  											alpha_i[alpha_index_i[active_size_i[i]]], minG))
  							{
  								swap(alpha_index_i[m], alpha_index_i[active_size_i[i]]);
  								swap(G[m], G[active_size_i[i]]);
  								if(y_index[i] == active_size_i[i])
  									y_index[i] = m;
  								else if(y_index[i] == m)
  									y_index[i] = active_size_i[i];
  								break;
  							}
  							active_size_i[i]--;
  						}
  					}
  				}
  
  				if(active_size_i[i] <= 1)
  				{
  					active_size--;
  					swap(index[s], index[active_size]);
  					s--;
  					continue;
  				}
  
  				if(maxG-minG <= 1e-12)
  					continue;
  				else
  					stopping = max(maxG - minG, stopping);
  
  				for(m=0;m<active_size_i[i];m++)
  					B[m] = G[m] - Ai*alpha_i[alpha_index_i[m]] ;
  
  				solve_sub_problem(Ai, y_index[i], C[GETI(i)], active_size_i[i], alpha_new);
  				int nz_d = 0;
  				for(m=0;m<active_size_i[i];m++)
  				{
  					double d = alpha_new[m] - alpha_i[alpha_index_i[m]];
  					alpha_i[alpha_index_i[m]] = alpha_new[m];
  					if(fabs(d) >= 1e-12)
  					{
  						d_ind[nz_d] = alpha_index_i[m];
  						d_val[nz_d] = d;
  						nz_d++;
  					}
  				}
  
  				xi = prob->x[i];
  				while(xi->index != -1)
  				{
  					double *w_i = &w[(xi->index-1)*nr_class];
  					for(m=0;m<nz_d;m++)
  						w_i[d_ind[m]] += d_val[m]*xi->value;
  					xi++;
  				}
  			}
  		}
  
  		iter++;
  		if(iter % 10 == 0)
  		{
  			info(".");
  		}
  
  		if(stopping < eps_shrink)
  		{
  			if(stopping < eps && start_from_all == true)
  				break;
  			else
  			{
  				active_size = l;
  				for(i=0;i<l;i++)
  					active_size_i[i] = nr_class;
  				info("*");
  				eps_shrink = max(eps_shrink/2, eps);
  				start_from_all = true;
  			}
  		}
  		else
  			start_from_all = false;
  	}
  
  	info("\noptimization finished, #iter = %d\n",iter);
  	if (iter >= max_iter)
  		info("\nWARNING: reaching max number of iterations\n");
  
  	// calculate objective value
  	double v = 0;
  	int nSV = 0;
  	for(i=0;i<w_size*nr_class;i++)
  		v += w[i]*w[i];
  	v = 0.5*v;
  	for(i=0;i<l*nr_class;i++)
  	{
  		v += alpha[i];
  		if(fabs(alpha[i]) > 0)
  			nSV++;
  	}
  	for(i=0;i<l;i++)
  		v -= alpha[i*nr_class+(int)prob->y[i]];
  	info("Objective value = %lf\n",v);
  	info("nSV = %d\n",nSV);
  
  	delete [] alpha;
  	delete [] alpha_new;
  	delete [] index;
  	delete [] QD;
  	delete [] d_ind;
  	delete [] d_val;
  	delete [] alpha_index;
  	delete [] y_index;
  	delete [] active_size_i;
  }
  
  // A coordinate descent algorithm for 
  // L1-loss and L2-loss SVM dual problems
  //
  //  min_\alpha  0.5(\alpha^T (Q + D)\alpha) - e^T \alpha,
  //    s.t.      0 <= \alpha_i <= upper_bound_i,
  // 
  //  where Qij = yi yj xi^T xj and
  //  D is a diagonal matrix 
  //
  // In L1-SVM case:
  // 		upper_bound_i = Cp if y_i = 1
  // 		upper_bound_i = Cn if y_i = -1
  // 		D_ii = 0
  // In L2-SVM case:
  // 		upper_bound_i = INF
  // 		D_ii = 1/(2*Cp)	if y_i = 1
  // 		D_ii = 1/(2*Cn)	if y_i = -1
  //
  // Given: 
  // x, y, Cp, Cn
  // eps is the stopping tolerance
  //
  // solution will be put in w
  // 
  // See Algorithm 3 of Hsieh et al., ICML 2008
  
  #undef GETI
  #define GETI(i) (y[i]+1)
  // To support weights for instances, use GETI(i) (i)
  
  static void solve_l2r_l1l2_svc(
  	const problem *prob, double *w, double eps,
  	double Cp, double Cn, int solver_type)
  {
  	int l = prob->l;
  	int w_size = prob->n;
  	int i, s, iter = 0;
  	double C, d, G;
  	double *QD = new double[l];
  	int max_iter = 1000;
  	int *index = new int[l];
  	double *alpha = new double[l];
  	schar *y = new schar[l];
  	int active_size = l;
  
  	// PG: projected gradient, for shrinking and stopping
  	double PG;
  	double PGmax_old = INF;
  	double PGmin_old = -INF;
  	double PGmax_new, PGmin_new;
  
  	// default solver_type: L2R_L2LOSS_SVC_DUAL
  	double diag[3] = {0.5/Cn, 0, 0.5/Cp};
  	double upper_bound[3] = {INF, 0, INF};
  	if(solver_type == L2R_L1LOSS_SVC_DUAL)
  	{
  		diag[0] = 0;
  		diag[2] = 0;
  		upper_bound[0] = Cn;
  		upper_bound[2] = Cp;
  	}
  
  	for(i=0; i<l; i++)
  	{
  		if(prob->y[i] > 0)
  		{
  			y[i] = +1;
  		}
  		else
  		{
  			y[i] = -1;
  		}
  	}
  
  	// Initial alpha can be set here. Note that
  	// 0 <= alpha[i] <= upper_bound[GETI(i)]
  	for(i=0; i<l; i++)
  		alpha[i] = 0;
  
  	for(i=0; i<w_size; i++)
  		w[i] = 0;
  	for(i=0; i<l; i++)
  	{
  		QD[i] = diag[GETI(i)];
  
  		feature_node *xi = prob->x[i];
  		while (xi->index != -1)
  		{
  			double val = xi->value;
  			QD[i] += val*val;
  			w[xi->index-1] += y[i]*alpha[i]*val;
  			xi++;
  		}
  		index[i] = i;
  	}
  
  	while (iter < max_iter)
  	{
  		PGmax_new = -INF;
  		PGmin_new = INF;
  
  		for (i=0; i<active_size; i++)
  		{
  			int j = i+rand()%(active_size-i);
  			swap(index[i], index[j]);
  		}
  
  		for (s=0; s<active_size; s++)
  		{
  			i = index[s];
  			G = 0;
  			schar yi = y[i];
  
  			feature_node *xi = prob->x[i];
  			while(xi->index!= -1)
  			{
  				G += w[xi->index-1]*(xi->value);
  				xi++;
  			}
  			G = G*yi-1;
  
  			C = upper_bound[GETI(i)];
  			G += alpha[i]*diag[GETI(i)];
  
  			PG = 0;
  			if (alpha[i] == 0)
  			{
  				if (G > PGmax_old)
  				{
  					active_size--;
  					swap(index[s], index[active_size]);
  					s--;
  					continue;
  				}
  				else if (G < 0)
  					PG = G;
  			}
  			else if (alpha[i] == C)
  			{
  				if (G < PGmin_old)
  				{
  					active_size--;
  					swap(index[s], index[active_size]);
  					s--;
  					continue;
  				}
  				else if (G > 0)
  					PG = G;
  			}
  			else
  				PG = G;
  
  			PGmax_new = max(PGmax_new, PG);
  			PGmin_new = min(PGmin_new, PG);
  
  			if(fabs(PG) > 1.0e-12)
  			{
  				double alpha_old = alpha[i];
  				alpha[i] = min(max(alpha[i] - G/QD[i], 0.0), C);
  				d = (alpha[i] - alpha_old)*yi;
  				xi = prob->x[i];
  				while (xi->index != -1)
  				{
  					w[xi->index-1] += d*xi->value;
  					xi++;
  				}
  			}
  		}
  
  		iter++;
  		if(iter % 10 == 0)
  			info(".");
  
  		if(PGmax_new - PGmin_new <= eps)
  		{
  			if(active_size == l)
  				break;
  			else
  			{
  				active_size = l;
  				info("*");
  				PGmax_old = INF;
  				PGmin_old = -INF;
  				continue;
  			}
  		}
  		PGmax_old = PGmax_new;
  		PGmin_old = PGmin_new;
  		if (PGmax_old <= 0)
  			PGmax_old = INF;
  		if (PGmin_old >= 0)
  			PGmin_old = -INF;
  	}
  
  	info("\noptimization finished, #iter = %d\n",iter);
  	if (iter >= max_iter)
  		info("\nWARNING: reaching max number of iterations\nUsing -s 2 may be faster (also see FAQ)\n\n");
  
  	// calculate objective value
  
  	double v = 0;
  	int nSV = 0;
  	for(i=0; i<w_size; i++)
  		v += w[i]*w[i];
  	for(i=0; i<l; i++)
  	{
  		v += alpha[i]*(alpha[i]*diag[GETI(i)] - 2);
  		if(alpha[i] > 0)
  			++nSV;
  	}
  	info("Objective value = %lf\n",v/2);
  	info("nSV = %d\n",nSV);
  
  	delete [] QD;
  	delete [] alpha;
  	delete [] y;
  	delete [] index;
  }
  
  
  // A coordinate descent algorithm for 
  // L1-loss and L2-loss epsilon-SVR dual problem
  //
  //  min_\beta  0.5\beta^T (Q + diag(lambda)) \beta - p \sum_{i=1}^l|\beta_i| + \sum_{i=1}^l yi\beta_i,
  //    s.t.      -upper_bound_i <= \beta_i <= upper_bound_i,
  // 
  //  where Qij = xi^T xj and
  //  D is a diagonal matrix 
  //
  // In L1-SVM case:
  // 		upper_bound_i = C
  // 		lambda_i = 0
  // In L2-SVM case:
  // 		upper_bound_i = INF
  // 		lambda_i = 1/(2*C)
  //
  // Given: 
  // x, y, p, C
  // eps is the stopping tolerance
  //
  // solution will be put in w
  //
  // See Algorithm 4 of Ho and Lin, 2012   
  
  #undef GETI
  #define GETI(i) (0)
  // To support weights for instances, use GETI(i) (i)
  
  static void solve_l2r_l1l2_svr(
  	const problem *prob, double *w, const parameter *param,
  	int solver_type)
  {
  	int l = prob->l;
  	double C = param->C;
  	double p = param->p;
  	int w_size = prob->n;
  	double eps = param->eps;
  	int i, s, iter = 0;
  	int max_iter = 1000;
  	int active_size = l;
  	int *index = new int[l];
  
  	double d, G, H;
  	double Gmax_old = INF;
  	double Gmax_new, Gnorm1_new;
  	double Gnorm1_init;
  	double *beta = new double[l];
  	double *QD = new double[l];
  	double *y = prob->y;
  
  	// L2R_L2LOSS_SVR_DUAL
  	double lambda[1], upper_bound[1];
  	lambda[0] = 0.5/C;
  	upper_bound[0] = INF;
  
  	if(solver_type == L2R_L1LOSS_SVR_DUAL)
  	{
  		lambda[0] = 0;
  		upper_bound[0] = C;
  	}
  
  	// Initial beta can be set here. Note that
  	// -upper_bound <= beta[i] <= upper_bound
  	for(i=0; i<l; i++)
  		beta[i] = 0;
  
  	for(i=0; i<w_size; i++)
  		w[i] = 0;
  	for(i=0; i<l; i++)
  	{
  		QD[i] = 0;
  		feature_node *xi = prob->x[i];
  		while(xi->index != -1)
  		{
  			double val = xi->value;
  			QD[i] += val*val;
  			w[xi->index-1] += beta[i]*val;
  			xi++;
  		}
  
  		index[i] = i;
  	}
  
  
  	while(iter < max_iter)
  	{
  		Gmax_new = 0;
  		Gnorm1_new = 0;
  
  		for(i=0; i<active_size; i++)
  		{
  			int j = i+rand()%(active_size-i);
  			swap(index[i], index[j]);
  		}
  
  		for(s=0; s<active_size; s++)
  		{
  			i = index[s];
  			G = -y[i] + lambda[GETI(i)]*beta[i];
  			H = QD[i] + lambda[GETI(i)];
  
  			feature_node *xi = prob->x[i];
  			while(xi->index != -1)
  			{
  				int ind = xi->index-1;
  				double val = xi->value;
  				G += val*w[ind];
  				xi++;
  			}
  
  			double Gp = G+p;
  			double Gn = G-p;
  			double violation = 0;
  			if(beta[i] == 0)
  			{
  				if(Gp < 0)
  					violation = -Gp;
  				else if(Gn > 0)
  					violation = Gn;
  				else if(Gp>Gmax_old && Gn<-Gmax_old)
  				{
  					active_size--;
  					swap(index[s], index[active_size]);
  					s--;
  					continue;
  				}
  			}
  			else if(beta[i] >= upper_bound[GETI(i)])
  			{
  				if(Gp > 0)
  					violation = Gp;
  				else if(Gp < -Gmax_old)
  				{
  					active_size--;
  					swap(index[s], index[active_size]);
  					s--;
  					continue;
  				}
  			}
  			else if(beta[i] <= -upper_bound[GETI(i)])
  			{
  				if(Gn < 0)
  					violation = -Gn;
  				else if(Gn > Gmax_old)
  				{
  					active_size--;
  					swap(index[s], index[active_size]);
  					s--;
  					continue;
  				}
  			}
  			else if(beta[i] > 0)
  				violation = fabs(Gp);
  			else
  				violation = fabs(Gn);
  
  			Gmax_new = max(Gmax_new, violation);
  			Gnorm1_new += violation;
  
  			// obtain Newton direction d
  			if(Gp < H*beta[i])
  				d = -Gp/H;
  			else if(Gn > H*beta[i])
  				d = -Gn/H;
  			else
  				d = -beta[i];
  
  			if(fabs(d) < 1.0e-12)
  				continue;
  
  			double beta_old = beta[i];
  			beta[i] = min(max(beta[i]+d, -upper_bound[GETI(i)]), upper_bound[GETI(i)]);
  			d = beta[i]-beta_old;
  
  			if(d != 0)
  			{
  				xi = prob->x[i];
  				while(xi->index != -1)
  				{
  					w[xi->index-1] += d*xi->value;
  					xi++;
  				}
  			}
  		}
  
  		if(iter == 0)
  			Gnorm1_init = Gnorm1_new;
  		iter++;
  		if(iter % 10 == 0)
  			info(".");
  
  		if(Gnorm1_new <= eps*Gnorm1_init)
  		{
  			if(active_size == l)
  				break;
  			else
  			{
  				active_size = l;
  				info("*");
  				Gmax_old = INF;
  				continue;
  			}
  		}
  
  		Gmax_old = Gmax_new;
  	}
  
  	info("\noptimization finished, #iter = %d\n", iter);
  	if(iter >= max_iter)
  		info("\nWARNING: reaching max number of iterations\nUsing -s 11 may be faster\n\n");
  
  	// calculate objective value
  	double v = 0;
  	int nSV = 0;
  	for(i=0; i<w_size; i++)
  		v += w[i]*w[i];
  	v = 0.5*v;
  	for(i=0; i<l; i++)
  	{
  		v += p*fabs(beta[i]) - y[i]*beta[i] + 0.5*lambda[GETI(i)]*beta[i]*beta[i];
  		if(beta[i] != 0)
  			nSV++;
  	}
  
  	info("Objective value = %lf\n", v);
  	info("nSV = %d\n",nSV);
  
  	delete [] beta;
  	delete [] QD;
  	delete [] index;
  }
  
  
  // A coordinate descent algorithm for 
  // the dual of L2-regularized logistic regression problems
  //
  //  min_\alpha  0.5(\alpha^T Q \alpha) + \sum \alpha_i log (\alpha_i) + (upper_bound_i - \alpha_i) log (upper_bound_i - \alpha_i),
  //    s.t.      0 <= \alpha_i <= upper_bound_i,
  // 
  //  where Qij = yi yj xi^T xj and 
  //  upper_bound_i = Cp if y_i = 1
  //  upper_bound_i = Cn if y_i = -1
  //
  // Given: 
  // x, y, Cp, Cn
  // eps is the stopping tolerance
  //
  // solution will be put in w
  //
  // See Algorithm 5 of Yu et al., MLJ 2010
  
  #undef GETI
  #define GETI(i) (y[i]+1)
  // To support weights for instances, use GETI(i) (i)
  
  void solve_l2r_lr_dual(const problem *prob, double *w, double eps, double Cp, double Cn)
  {
  	int l = prob->l;
  	int w_size = prob->n;
  	int i, s, iter = 0;
  	double *xTx = new double[l];
  	int max_iter = 1000;
  	int *index = new int[l];	
  	double *alpha = new double[2*l]; // store alpha and C - alpha
  	schar *y = new schar[l];
  	int max_inner_iter = 100; // for inner Newton
  	double innereps = 1e-2;
  	double innereps_min = min(1e-8, eps);
  	double upper_bound[3] = {Cn, 0, Cp};
  
  	for(i=0; i<l; i++)
  	{
  		if(prob->y[i] > 0)
  		{
  			y[i] = +1;
  		}
  		else
  		{
  			y[i] = -1;
  		}
  	}
  	
  	// Initial alpha can be set here. Note that
  	// 0 < alpha[i] < upper_bound[GETI(i)]
  	// alpha[2*i] + alpha[2*i+1] = upper_bound[GETI(i)]
  	for(i=0; i<l; i++)
  	{
  		alpha[2*i] = min(0.001*upper_bound[GETI(i)], 1e-8);
  		alpha[2*i+1] = upper_bound[GETI(i)] - alpha[2*i];
  	}
  
  	for(i=0; i<w_size; i++)
  		w[i] = 0;
  	for(i=0; i<l; i++)
  	{
  		xTx[i] = 0;
  		feature_node *xi = prob->x[i];
  		while (xi->index != -1)
  		{
  			double val = xi->value;
  			xTx[i] += val*val;
  			w[xi->index-1] += y[i]*alpha[2*i]*val;
  			xi++;
  		}
  		index[i] = i;
  	}
  
  	while (iter < max_iter)
  	{
  		for (i=0; i<l; i++)
  		{
  			int j = i+rand()%(l-i);
  			swap(index[i], index[j]);
  		}
  		int newton_iter = 0;
  		double Gmax = 0;
  		for (s=0; s<l; s++)
  		{
  			i = index[s];
  			schar yi = y[i];
  			double C = upper_bound[GETI(i)];
  			double ywTx = 0, xisq = xTx[i];
  			feature_node *xi = prob->x[i];
  			while (xi->index != -1)
  			{
  				ywTx += w[xi->index-1]*xi->value;
  				xi++;
  			}
  			ywTx *= y[i];
  			double a = xisq, b = ywTx;
  
  			// Decide to minimize g_1(z) or g_2(z)
  			int ind1 = 2*i, ind2 = 2*i+1, sign = 1;
  			if(0.5*a*(alpha[ind2]-alpha[ind1])+b < 0)
  			{
  				ind1 = 2*i+1;
  				ind2 = 2*i;
  				sign = -1;
  			}
  
  			//  g_t(z) = z*log(z) + (C-z)*log(C-z) + 0.5a(z-alpha_old)^2 + sign*b(z-alpha_old)
  			double alpha_old = alpha[ind1];
  			double z = alpha_old;
  			if(C - z < 0.5 * C)
  				z = 0.1*z;
  			double gp = a*(z-alpha_old)+sign*b+log(z/(C-z));
  			Gmax = max(Gmax, fabs(gp));
  
  			// Newton method on the sub-problem
  			const double eta = 0.1; // xi in the paper
  			int inner_iter = 0;
  			while (inner_iter <= max_inner_iter)
  			{
  				if(fabs(gp) < innereps)
  					break;
  				double gpp = a + C/(C-z)/z;
  				double tmpz = z - gp/gpp;
  				if(tmpz <= 0)
  					z *= eta;
  				else // tmpz in (0, C)
  					z = tmpz;
  				gp = a*(z-alpha_old)+sign*b+log(z/(C-z));
  				newton_iter++;
  				inner_iter++;
  			}
  
  			if(inner_iter > 0) // update w
  			{
  				alpha[ind1] = z;
  				alpha[ind2] = C-z;
  				xi = prob->x[i];
  				while (xi->index != -1)
  				{
  					w[xi->index-1] += sign*(z-alpha_old)*yi*xi->value;
  					xi++;
  				}
  			}
  		}
  
  		iter++;
  		if(iter % 10 == 0)
  			info(".");
  
  		if(Gmax < eps)
  			break;
  
  		if(newton_iter <= l/10)
  			innereps = max(innereps_min, 0.1*innereps);
  
  	}
  
  	info("\noptimization finished, #iter = %d\n",iter);
  	if (iter >= max_iter)
  		info("\nWARNING: reaching max number of iterations\nUsing -s 0 may be faster (also see FAQ)\n\n");
  
  	// calculate objective value
  
  	double v = 0;
  	for(i=0; i<w_size; i++)
  		v += w[i] * w[i];
  	v *= 0.5;
  	for(i=0; i<l; i++)
  		v += alpha[2*i] * log(alpha[2*i]) + alpha[2*i+1] * log(alpha[2*i+1])
  			- upper_bound[GETI(i)] * log(upper_bound[GETI(i)]);
  	info("Objective value = %lf\n", v);
  
  	delete [] xTx;
  	delete [] alpha;
  	delete [] y;
  	delete [] index;
  }
  
  // A coordinate descent algorithm for 
  // L1-regularized L2-loss support vector classification
  //
  //  min_w \sum |wj| + C \sum max(0, 1-yi w^T xi)^2,
  //
  // Given: 
  // x, y, Cp, Cn
  // eps is the stopping tolerance
  //
  // solution will be put in w
  //
  // See Yuan et al. (2010) and appendix of LIBLINEAR paper, Fan et al. (2008)
  
  #undef GETI
  #define GETI(i) (y[i]+1)
  // To support weights for instances, use GETI(i) (i)
  
  static void solve_l1r_l2_svc(
  	problem *prob_col, double *w, double eps,
  	double Cp, double Cn)
  {
  	int l = prob_col->l;
  	int w_size = prob_col->n;
  	int j, s, iter = 0;
  	int max_iter = 1000;
  	int active_size = w_size;
  	int max_num_linesearch = 20;
  
  	double sigma = 0.01;
  	double d, G_loss, G, H;
  	double Gmax_old = INF;
  	double Gmax_new, Gnorm1_new;
  	double Gnorm1_init;
  	double d_old, d_diff;
  	double loss_old, loss_new;
  	double appxcond, cond;
  
  	int *index = new int[w_size];
  	schar *y = new schar[l];
  	double *b = new double[l]; // b = 1-ywTx
  	double *xj_sq = new double[w_size];
  	feature_node *x;
  
  	double C[3] = {Cn,0,Cp};
  
  	// Initial w can be set here.
  	for(j=0; j<w_size; j++)
  		w[j] = 0;
  
  	for(j=0; j<l; j++)
  	{
  		b[j] = 1;
  		if(prob_col->y[j] > 0)
  			y[j] = 1;
  		else
  			y[j] = -1;
  	}
  	for(j=0; j<w_size; j++)
  	{
  		index[j] = j;
  		xj_sq[j] = 0;
  		x = prob_col->x[j];
  		while(x->index != -1)
  		{
  			int ind = x->index-1;
  			x->value *= y[ind]; // x->value stores yi*xij
  			double val = x->value;
  			b[ind] -= w[j]*val;
  			xj_sq[j] += C[GETI(ind)]*val*val;
  			x++;
  		}
  	}
  
  	while(iter < max_iter)
  	{
  		Gmax_new = 0;
  		Gnorm1_new = 0;
  
  		for(j=0; j<active_size; j++)
  		{
  			int i = j+rand()%(active_size-j);
  			swap(index[i], index[j]);
  		}
  
  		for(s=0; s<active_size; s++)
  		{
  			j = index[s];
  			G_loss = 0;
  			H = 0;
  
  			x = prob_col->x[j];
  			while(x->index != -1)
  			{
  				int ind = x->index-1;
  				if(b[ind] > 0)
  				{
  					double val = x->value;
  					double tmp = C[GETI(ind)]*val;
  					G_loss -= tmp*b[ind];
  					H += tmp*val;
  				}
  				x++;
  			}
  			G_loss *= 2;
  
  			G = G_loss;
  			H *= 2;
  			H = max(H, 1e-12);
  
  			double Gp = G+1;
  			double Gn = G-1;
  			double violation = 0;
  			if(w[j] == 0)
  			{
  				if(Gp < 0)
  					violation = -Gp;
  				else if(Gn > 0)
  					violation = Gn;
  				else if(Gp>Gmax_old/l && Gn<-Gmax_old/l)
  				{
  					active_size--;
  					swap(index[s], index[active_size]);
  					s--;
  					continue;
  				}
  			}
  			else if(w[j] > 0)
  				violation = fabs(Gp);
  			else
  				violation = fabs(Gn);
  
  			Gmax_new = max(Gmax_new, violation);
  			Gnorm1_new += violation;
  
  			// obtain Newton direction d
  			if(Gp < H*w[j])
  				d = -Gp/H;
  			else if(Gn > H*w[j])
  				d = -Gn/H;
  			else
  				d = -w[j];
  
  			if(fabs(d) < 1.0e-12)
  				continue;
  
  			double delta = fabs(w[j]+d)-fabs(w[j]) + G*d;
  			d_old = 0;
  			int num_linesearch;
  			for(num_linesearch=0; num_linesearch < max_num_linesearch; num_linesearch++)
  			{
  				d_diff = d_old - d;
  				cond = fabs(w[j]+d)-fabs(w[j]) - sigma*delta;
  
  				appxcond = xj_sq[j]*d*d + G_loss*d + cond;
  				if(appxcond <= 0)
  				{
  					x = prob_col->x[j];
  					while(x->index != -1)
  					{
  						b[x->index-1] += d_diff*x->value;
  						x++;
  					}
  					break;
  				}
  
  				if(num_linesearch == 0)
  				{
  					loss_old = 0;
  					loss_new = 0;
  					x = prob_col->x[j];
  					while(x->index != -1)
  					{
  						int ind = x->index-1;
  						if(b[ind] > 0)
  							loss_old += C[GETI(ind)]*b[ind]*b[ind];
  						double b_new = b[ind] + d_diff*x->value;
  						b[ind] = b_new;
  						if(b_new > 0)
  							loss_new += C[GETI(ind)]*b_new*b_new;
  						x++;
  					}
  				}
  				else
  				{
  					loss_new = 0;
  					x = prob_col->x[j];
  					while(x->index != -1)
  					{
  						int ind = x->index-1;
  						double b_new = b[ind] + d_diff*x->value;
  						b[ind] = b_new;
  						if(b_new > 0)
  							loss_new += C[GETI(ind)]*b_new*b_new;
  						x++;
  					}
  				}
  
  				cond = cond + loss_new - loss_old;
  				if(cond <= 0)
  					break;
  				else
  				{
  					d_old = d;
  					d *= 0.5;
  					delta *= 0.5;
  				}
  			}
  
  			w[j] += d;
  
  			// recompute b[] if line search takes too many steps
  			if(num_linesearch >= max_num_linesearch)
  			{
  				info("#");
  				for(int i=0; i<l; i++)
  					b[i] = 1;
  
  				for(int i=0; i<w_size; i++)
  				{
  					if(w[i]==0) continue;
  					x = prob_col->x[i];
  					while(x->index != -1)
  					{
  						b[x->index-1] -= w[i]*x->value;
  						x++;
  					}
  				}
  			}
  		}
  
  		if(iter == 0)
  			Gnorm1_init = Gnorm1_new;
  		iter++;
  		if(iter % 10 == 0)
  			info(".");
  
  		if(Gnorm1_new <= eps*Gnorm1_init)
  		{
  			if(active_size == w_size)
  				break;
  			else
  			{
  				active_size = w_size;
  				info("*");
  				Gmax_old = INF;
  				continue;
  			}
  		}
  
  		Gmax_old = Gmax_new;
  	}
  
  	info("\noptimization finished, #iter = %d\n", iter);
  	if(iter >= max_iter)
  		info("\nWARNING: reaching max number of iterations\n");
  
  	// calculate objective value
  
  	double v = 0;
  	int nnz = 0;
  	for(j=0; j<w_size; j++)
  	{
  		x = prob_col->x[j];
  		while(x->index != -1)
  		{
  			x->value *= prob_col->y[x->index-1]; // restore x->value
  			x++;
  		}
  		if(w[j] != 0)
  		{
  			v += fabs(w[j]);
  			nnz++;
  		}
  	}
  	for(j=0; j<l; j++)
  		if(b[j] > 0)
  			v += C[GETI(j)]*b[j]*b[j];
  
  	info("Objective value = %lf\n", v);
  	info("#nonzeros/#features = %d/%d\n", nnz, w_size);
  
  	delete [] index;
  	delete [] y;
  	delete [] b;
  	delete [] xj_sq;
  }
  
  // A coordinate descent algorithm for 
  // L1-regularized logistic regression problems
  //
  //  min_w \sum |wj| + C \sum log(1+exp(-yi w^T xi)),
  //
  // Given: 
  // x, y, Cp, Cn
  // eps is the stopping tolerance
  //
  // solution will be put in w
  //
  // See Yuan et al. (2011) and appendix of LIBLINEAR paper, Fan et al. (2008)
  
  #undef GETI
  #define GETI(i) (y[i]+1)
  // To support weights for instances, use GETI(i) (i)
  
  static void solve_l1r_lr(
  	const problem *prob_col, double *w, double eps,
  	double Cp, double Cn)
  {
  	int l = prob_col->l;
  	int w_size = prob_col->n;
  	int j, s, newton_iter=0, iter=0;
  	int max_newton_iter = 100;
  	int max_iter = 1000;
  	int max_num_linesearch = 20;
  	int active_size;
  	int QP_active_size;
  
  	double nu = 1e-12;
  	double inner_eps = 1;
  	double sigma = 0.01;
  	double w_norm, w_norm_new;
  	double z, G, H;
  	double Gnorm1_init;
  	double Gmax_old = INF;
  	double Gmax_new, Gnorm1_new;
  	double QP_Gmax_old = INF;
  	double QP_Gmax_new, QP_Gnorm1_new;
  	double delta, negsum_xTd, cond;
  
  	int *index = new int[w_size];
  	schar *y = new schar[l];
  	double *Hdiag = new double[w_size];
  	double *Grad = new double[w_size];
  	double *wpd = new double[w_size];
  	double *xjneg_sum = new double[w_size];
  	double *xTd = new double[l];
  	double *exp_wTx = new double[l];
  	double *exp_wTx_new = new double[l];
  	double *tau = new double[l];
  	double *D = new double[l];
  	feature_node *x;
  
  	double C[3] = {Cn,0,Cp};
  
  	// Initial w can be set here.
  	for(j=0; j<w_size; j++)
  		w[j] = 0;
  
  	for(j=0; j<l; j++)
  	{
  		if(prob_col->y[j] > 0)
  			y[j] = 1;
  		else
  			y[j] = -1;
  
  		exp_wTx[j] = 0;
  	}
  
  	w_norm = 0;
  	for(j=0; j<w_size; j++)
  	{
  		w_norm += fabs(w[j]);
  		wpd[j] = w[j];
  		index[j] = j;
  		xjneg_sum[j] = 0;
  		x = prob_col->x[j];
  		while(x->index != -1)
  		{
  			int ind = x->index-1;
  			double val = x->value;
  			exp_wTx[ind] += w[j]*val;
  			if(y[ind] == -1)
  				xjneg_sum[j] += C[GETI(ind)]*val;
  			x++;
  		}
  	}
  	for(j=0; j<l; j++)
  	{
  		exp_wTx[j] = exp(exp_wTx[j]);
  		double tau_tmp = 1/(1+exp_wTx[j]);
  		tau[j] = C[GETI(j)]*tau_tmp;
  		D[j] = C[GETI(j)]*exp_wTx[j]*tau_tmp*tau_tmp;
  	}
  
  	while(newton_iter < max_newton_iter)
  	{
  		Gmax_new = 0;
  		Gnorm1_new = 0;
  		active_size = w_size;
  
  		for(s=0; s<active_size; s++)
  		{
  			j = index[s];
  			Hdiag[j] = nu;
  			Grad[j] = 0;
  
  			double tmp = 0;
  			x = prob_col->x[j];
  			while(x->index != -1)
  			{
  				int ind = x->index-1;
  				Hdiag[j] += x->value*x->value*D[ind];
  				tmp += x->value*tau[ind];
  				x++;
  			}
  			Grad[j] = -tmp + xjneg_sum[j];
  
  			double Gp = Grad[j]+1;
  			double Gn = Grad[j]-1;
  			double violation = 0;
  			if(w[j] == 0)
  			{
  				if(Gp < 0)
  					violation = -Gp;
  				else if(Gn > 0)
  					violation = Gn;
  				//outer-level shrinking
  				else if(Gp>Gmax_old/l && Gn<-Gmax_old/l)
  				{
  					active_size--;
  					swap(index[s], index[active_size]);
  					s--;
  					continue;
  				}
  			}
  			else if(w[j] > 0)
  				violation = fabs(Gp);
  			else
  				violation = fabs(Gn);
  
  			Gmax_new = max(Gmax_new, violation);
  			Gnorm1_new += violation;
  		}
  
  		if(newton_iter == 0)
  			Gnorm1_init = Gnorm1_new;
  
  		if(Gnorm1_new <= eps*Gnorm1_init)
  			break;
  
  		iter = 0;
  		QP_Gmax_old = INF;
  		QP_active_size = active_size;
  
  		for(int i=0; i<l; i++)
  			xTd[i] = 0;
  
  		// optimize QP over wpd
  		while(iter < max_iter)
  		{
  			QP_Gmax_new = 0;
  			QP_Gnorm1_new = 0;
  
  			for(j=0; j<QP_active_size; j++)
  			{
  				int i = j+rand()%(QP_active_size-j);
  				swap(index[i], index[j]);
  			}
  
  			for(s=0; s<QP_active_size; s++)
  			{
  				j = index[s];
  				H = Hdiag[j];
  
  				x = prob_col->x[j];
  				G = Grad[j] + (wpd[j]-w[j])*nu;
  				while(x->index != -1)
  				{
  					int ind = x->index-1;
  					G += x->value*D[ind]*xTd[ind];
  					x++;
  				}
  
  				double Gp = G+1;
  				double Gn = G-1;
  				double violation = 0;
  				if(wpd[j] == 0)
  				{
  					if(Gp < 0)
  						violation = -Gp;
  					else if(Gn > 0)
  						violation = Gn;
  					//inner-level shrinking
  					else if(Gp>QP_Gmax_old/l && Gn<-QP_Gmax_old/l)
  					{
  						QP_active_size--;
  						swap(index[s], index[QP_active_size]);
  						s--;
  						continue;
  					}
  				}
  				else if(wpd[j] > 0)
  					violation = fabs(Gp);
  				else
  					violation = fabs(Gn);
  
  				QP_Gmax_new = max(QP_Gmax_new, violation);
  				QP_Gnorm1_new += violation;
  
  				// obtain solution of one-variable problem
  				if(Gp < H*wpd[j])
  					z = -Gp/H;
  				else if(Gn > H*wpd[j])
  					z = -Gn/H;
  				else
  					z = -wpd[j];
  
  				if(fabs(z) < 1.0e-12)
  					continue;
  				z = min(max(z,-10.0),10.0);
  
  				wpd[j] += z;
  
  				x = prob_col->x[j];
  				while(x->index != -1)
  				{
  					int ind = x->index-1;
  					xTd[ind] += x->value*z;
  					x++;
  				}
  			}
  
  			iter++;
  
  			if(QP_Gnorm1_new <= inner_eps*Gnorm1_init)
  			{
  				//inner stopping
  				if(QP_active_size == active_size)
  					break;
  				//active set reactivation
  				else
  				{
  					QP_active_size = active_size;
  					QP_Gmax_old = INF;
  					continue;
  				}
  			}
  
  			QP_Gmax_old = QP_Gmax_new;
  		}
  
  		if(iter >= max_iter)
  			info("WARNING: reaching max number of inner iterations\n");
  
  		delta = 0;
  		w_norm_new = 0;
  		for(j=0; j<w_size; j++)
  		{
  			delta += Grad[j]*(wpd[j]-w[j]);
  			if(wpd[j] != 0)
  				w_norm_new += fabs(wpd[j]);
  		}
  		delta += (w_norm_new-w_norm);
  
  		negsum_xTd = 0;
  		for(int i=0; i<l; i++)
  			if(y[i] == -1)
  				negsum_xTd += C[GETI(i)]*xTd[i];
  
  		int num_linesearch;
  		for(num_linesearch=0; num_linesearch < max_num_linesearch; num_linesearch++)
  		{
  			cond = w_norm_new - w_norm + negsum_xTd - sigma*delta;
  
  			for(int i=0; i<l; i++)
  			{
  				double exp_xTd = exp(xTd[i]);
  				exp_wTx_new[i] = exp_wTx[i]*exp_xTd;
  				cond += C[GETI(i)]*log((1+exp_wTx_new[i])/(exp_xTd+exp_wTx_new[i]));
  			}
  
  			if(cond <= 0)
  			{
  				w_norm = w_norm_new;
  				for(j=0; j<w_size; j++)
  					w[j] = wpd[j];
  				for(int i=0; i<l; i++)
  				{
  					exp_wTx[i] = exp_wTx_new[i];
  					double tau_tmp = 1/(1+exp_wTx[i]);
  					tau[i] = C[GETI(i)]*tau_tmp;
  					D[i] = C[GETI(i)]*exp_wTx[i]*tau_tmp*tau_tmp;
  				}
  				break;
  			}
  			else
  			{
  				w_norm_new = 0;
  				for(j=0; j<w_size; j++)
  				{
  					wpd[j] = (w[j]+wpd[j])*0.5;
  					if(wpd[j] != 0)
  						w_norm_new += fabs(wpd[j]);
  				}
  				delta *= 0.5;
  				negsum_xTd *= 0.5;
  				for(int i=0; i<l; i++)
  					xTd[i] *= 0.5;
  			}
  		}
  
  		// Recompute some info due to too many line search steps
  		if(num_linesearch >= max_num_linesearch)
  		{
  			for(int i=0; i<l; i++)
  				exp_wTx[i] = 0;
  
  			for(int i=0; i<w_size; i++)
  			{
  				if(w[i]==0) continue;
  				x = prob_col->x[i];
  				while(x->index != -1)
  				{
  					exp_wTx[x->index-1] += w[i]*x->value;
  					x++;
  				}
  			}
  
  			for(int i=0; i<l; i++)
  				exp_wTx[i] = exp(exp_wTx[i]);
  		}
  
  		if(iter == 1)
  			inner_eps *= 0.25;
  
  		newton_iter++;
  		Gmax_old = Gmax_new;
  
  		info("iter %3d  #CD cycles %d\n", newton_iter, iter);
  	}
  
  	info("=========================\n");
  	info("optimization finished, #iter = %d\n", newton_iter);
  	if(newton_iter >= max_newton_iter)
  		info("WARNING: reaching max number of iterations\n");
  
  	// calculate objective value
  
  	double v = 0;
  	int nnz = 0;
  	for(j=0; j<w_size; j++)
  		if(w[j] != 0)
  		{
  			v += fabs(w[j]);
  			nnz++;
  		}
  	for(j=0; j<l; j++)
  		if(y[j] == 1)
  			v += C[GETI(j)]*log(1+1/exp_wTx[j]);
  		else
  			v += C[GETI(j)]*log(1+exp_wTx[j]);
  
  	info("Objective value = %lf\n", v);
  	info("#nonzeros/#features = %d/%d\n", nnz, w_size);
  
  	delete [] index;
  	delete [] y;
  	delete [] Hdiag;
  	delete [] Grad;
  	delete [] wpd;
  	delete [] xjneg_sum;
  	delete [] xTd;
  	delete [] exp_wTx;
  	delete [] exp_wTx_new;
  	delete [] tau;
  	delete [] D;
  }
  
  // transpose matrix X from row format to column format
  static void transpose(const problem *prob, feature_node **x_space_ret, problem *prob_col)
  {
  	int i;
  	int l = prob->l;
  	int n = prob->n;
  	long int nnz = 0;
  	long int *col_ptr = new long int [n+1];
  	feature_node *x_space;
  	prob_col->l = l;
  	prob_col->n = n;
  	prob_col->y = new double[l];
  	prob_col->x = new feature_node*[n];
  
  	for(i=0; i<l; i++)
  		prob_col->y[i] = prob->y[i];
  
  	for(i=0; i<n+1; i++)
  		col_ptr[i] = 0;
  	for(i=0; i<l; i++)
  	{
  		feature_node *x = prob->x[i];
  		while(x->index != -1)
  		{
  			nnz++;
  			col_ptr[x->index]++;
  			x++;
  		}
  	}
  	for(i=1; i<n+1; i++)
  		col_ptr[i] += col_ptr[i-1] + 1;
  
  	x_space = new feature_node[nnz+n];
  	for(i=0; i<n; i++)
  		prob_col->x[i] = &x_space[col_ptr[i]];
  
  	for(i=0; i<l; i++)
  	{
  		feature_node *x = prob->x[i];
  		while(x->index != -1)
  		{
  			int ind = x->index-1;
  			x_space[col_ptr[ind]].index = i+1; // starts from 1
  			x_space[col_ptr[ind]].value = x->value;
  			col_ptr[ind]++;
  			x++;
  		}
  	}
  	for(i=0; i<n; i++)
  		x_space[col_ptr[i]].index = -1;
  
  	*x_space_ret = x_space;
  
  	delete [] col_ptr;
  }
  
  // label: label name, start: begin of each class, count: #data of classes, perm: indices to the original data
  // perm, length l, must be allocated before calling this subroutine
  static void group_classes(const problem *prob, int *nr_class_ret, int **label_ret, int **start_ret, int **count_ret, int *perm)
  {
  	int l = prob->l;
  	int max_nr_class = 16;
  	int nr_class = 0;
  	int *label = Malloc(int,max_nr_class);
  	int *count = Malloc(int,max_nr_class);
  	int *data_label = Malloc(int,l);
  	int i;
  
  	for(i=0;i<l;i++)
  	{
  		int this_label = (int)prob->y[i];
  		int j;
  		for(j=0;j<nr_class;j++)
  		{
  			if(this_label == label[j])
  			{
  				++count[j];
  				break;
  			}
  		}
  		data_label[i] = j;
  		if(j == nr_class)
  		{
  			if(nr_class == max_nr_class)
  			{
  				max_nr_class *= 2;
  				label = (int *)realloc(label,max_nr_class*sizeof(int));
  				count = (int *)realloc(count,max_nr_class*sizeof(int));
  			}
  			label[nr_class] = this_label;
  			count[nr_class] = 1;
  			++nr_class;
  		}
  	}
  
  	//
  	// Labels are ordered by their first occurrence in the training set. 
  	// However, for two-class sets with -1/+1 labels and -1 appears first, 
  	// we swap labels to ensure that internally the binary SVM has positive data corresponding to the +1 instances.
  	//
  	if (nr_class == 2 && label[0] == -1 && label[1] == 1)
  	{
  		swap(label[0],label[1]);
  		swap(count[0],count[1]);
  		for(i=0;i<l;i++)
  		{
  			if(data_label[i] == 0)
  				data_label[i] = 1;
  			else
  				data_label[i] = 0;
  		}
  	}
  
  	int *start = Malloc(int,nr_class);
  	start[0] = 0;
  	for(i=1;i<nr_class;i++)
  		start[i] = start[i-1]+count[i-1];
  	for(i=0;i<l;i++)
  	{
  		perm[start[data_label[i]]] = i;
  		++start[data_label[i]];
  	}
  	start[0] = 0;
  	for(i=1;i<nr_class;i++)
  		start[i] = start[i-1]+count[i-1];
  
  	*nr_class_ret = nr_class;
  	*label_ret = label;
  	*start_ret = start;
  	*count_ret = count;
  	free(data_label);
  }
  
  static void train_one(const problem *prob, const parameter *param, double *w, double Cp, double Cn)
  {
  	double eps=param->eps;
  	int pos = 0;
  	int neg = 0;
  	for(int i=0;i<prob->l;i++)
  		if(prob->y[i] > 0)
  			pos++;
  	neg = prob->l - pos;
  
  	double primal_solver_tol = eps*max(min(pos,neg), 1)/prob->l;
  
  	function *fun_obj=NULL;
  	switch(param->solver_type)
  	{
  		case L2R_LR:
  		{
  			double *C = new double[prob->l];
  			for(int i = 0; i < prob->l; i++)
  			{
  				if(prob->y[i] > 0)
  					C[i] = Cp;
  				else
  					C[i] = Cn;
  			}
  			fun_obj=new l2r_lr_fun(prob, C);
  			TRON tron_obj(fun_obj, primal_solver_tol);
  			tron_obj.set_print_string(liblinear_print_string);
  			tron_obj.tron(w);
  			delete fun_obj;
  			delete[] C;
  			break;
  		}
  		case L2R_L2LOSS_SVC:
  		{
  			double *C = new double[prob->l];
  			for(int i = 0; i < prob->l; i++)
  			{
  				if(prob->y[i] > 0)
  					C[i] = Cp;
  				else
  					C[i] = Cn;
  			}
  			fun_obj=new l2r_l2_svc_fun(prob, C);
  			TRON tron_obj(fun_obj, primal_solver_tol);
  			tron_obj.set_print_string(liblinear_print_string);
  			tron_obj.tron(w);
  			delete fun_obj;
  			delete[] C;
  			break;
  		}
  		case L2R_L2LOSS_SVC_DUAL:
  			solve_l2r_l1l2_svc(prob, w, eps, Cp, Cn, L2R_L2LOSS_SVC_DUAL);
  			break;
  		case L2R_L1LOSS_SVC_DUAL:
  			solve_l2r_l1l2_svc(prob, w, eps, Cp, Cn, L2R_L1LOSS_SVC_DUAL);
  			break;
  		case L1R_L2LOSS_SVC:
  		{
  			problem prob_col;
  			feature_node *x_space = NULL;
  			transpose(prob, &x_space ,&prob_col);
  			solve_l1r_l2_svc(&prob_col, w, primal_solver_tol, Cp, Cn);
  			delete [] prob_col.y;
  			delete [] prob_col.x;
  			delete [] x_space;
  			break;
  		}
  		case L1R_LR:
  		{
  			problem prob_col;
  			feature_node *x_space = NULL;
  			transpose(prob, &x_space ,&prob_col);
  			solve_l1r_lr(&prob_col, w, primal_solver_tol, Cp, Cn);
  			delete [] prob_col.y;
  			delete [] prob_col.x;
  			delete [] x_space;
  			break;
  		}
  		case L2R_LR_DUAL:
  			solve_l2r_lr_dual(prob, w, eps, Cp, Cn);
  			break;
  		case L2R_L2LOSS_SVR:
  		{
  			double *C = new double[prob->l];
  			for(int i = 0; i < prob->l; i++)
  				C[i] = param->C;
  
  			fun_obj=new l2r_l2_svr_fun(prob, C, param->p);
  			TRON tron_obj(fun_obj, param->eps);
  			tron_obj.set_print_string(liblinear_print_string);
  			tron_obj.tron(w);
  			delete fun_obj;
  			delete[] C;
  			break;
  
  		}
  		case L2R_L1LOSS_SVR_DUAL:
  			solve_l2r_l1l2_svr(prob, w, param, L2R_L1LOSS_SVR_DUAL);
  			break;
  		case L2R_L2LOSS_SVR_DUAL:
  			solve_l2r_l1l2_svr(prob, w, param, L2R_L2LOSS_SVR_DUAL);
  			break;
  		default:
  			fprintf(stderr, "ERROR: unknown solver_type\n");
  			break;
  	}
  }
  
  //
  // Interface functions
  //
  model* train(const problem *prob, const parameter *param)
  {
  	int i,j;
  	int l = prob->l;
  	int n = prob->n;
  	int w_size = prob->n;
  	model *model_ = Malloc(model,1);
  
  	if(prob->bias>=0)
  		model_->nr_feature=n-1;
  	else
  		model_->nr_feature=n;
  	model_->param = *param;
  	model_->bias = prob->bias;
  
  	if(param->solver_type == L2R_L2LOSS_SVR ||
  	   param->solver_type == L2R_L1LOSS_SVR_DUAL ||
  	   param->solver_type == L2R_L2LOSS_SVR_DUAL)
  	{
  		model_->w = Malloc(double, w_size);
  		model_->nr_class = 2;
  		model_->label = NULL;
  		train_one(prob, param, &model_->w[0], 0, 0);
  	}
  	else
  	{
  		int nr_class;
  		int *label = NULL;
  		int *start = NULL;
  		int *count = NULL;
  		int *perm = Malloc(int,l);
  
  		// group training data of the same class
  		group_classes(prob,&nr_class,&label,&start,&count,perm);
  
  		model_->nr_class=nr_class;
  		model_->label = Malloc(int,nr_class);
  		for(i=0;i<nr_class;i++)
  			model_->label[i] = label[i];
  
  		// calculate weighted C
  		double *weighted_C = Malloc(double, nr_class);
  		for(i=0;i<nr_class;i++)
  			weighted_C[i] = param->C;
  		for(i=0;i<param->nr_weight;i++)
  		{
  			for(j=0;j<nr_class;j++)
  				if(param->weight_label[i] == label[j])
  					break;
  			if(j == nr_class)
  				fprintf(stderr,"WARNING: class label %d specified in weight is not found\n", param->weight_label[i]);
  			else
  				weighted_C[j] *= param->weight[i];
  		}
  
  		// constructing the subproblem
  		feature_node **x = Malloc(feature_node *,l);
  		for(i=0;i<l;i++)
  			x[i] = prob->x[perm[i]];
  
  		int k;
  		problem sub_prob;
  		sub_prob.l = l;
  		sub_prob.n = n;
  		sub_prob.x = Malloc(feature_node *,sub_prob.l);
  		sub_prob.y = Malloc(double,sub_prob.l);
  
  		for(k=0; k<sub_prob.l; k++)
  			sub_prob.x[k] = x[k];
  
  		// multi-class svm by Crammer and Singer
  		if(param->solver_type == MCSVM_CS)
  		{
  			model_->w=Malloc(double, n*nr_class);
  			for(i=0;i<nr_class;i++)
  				for(j=start[i];j<start[i]+count[i];j++)
  					sub_prob.y[j] = i;
  			Solver_MCSVM_CS Solver(&sub_prob, nr_class, weighted_C, param->eps);
  			Solver.Solve(model_->w);
  		}
  		else
  		{
  			if(nr_class == 2)
  			{
  				model_->w=Malloc(double, w_size);
  
  				int e0 = start[0]+count[0];
  				k=0;
  				for(; k<e0; k++)
  					sub_prob.y[k] = +1;
  				for(; k<sub_prob.l; k++)
  					sub_prob.y[k] = -1;
  
  				train_one(&sub_prob, param, &model_->w[0], weighted_C[0], weighted_C[1]);
  			}
  			else
  			{
  				model_->w=Malloc(double, w_size*nr_class);
  				double *w=Malloc(double, w_size);
  				for(i=0;i<nr_class;i++)
  				{
  					int si = start[i];
  					int ei = si+count[i];
  
  					k=0;
  					for(; k<si; k++)
  						sub_prob.y[k] = -1;
  					for(; k<ei; k++)
  						sub_prob.y[k] = +1;
  					for(; k<sub_prob.l; k++)
  						sub_prob.y[k] = -1;
  
  					train_one(&sub_prob, param, w, weighted_C[i], param->C);
  
  					for(int j=0;j<w_size;j++)
  						model_->w[j*nr_class+i] = w[j];
  				}
  				free(w);
  			}
  
  		}
  
  		free(x);
  		free(label);
  		free(start);
  		free(count);
  		free(perm);
  		free(sub_prob.x);
  		free(sub_prob.y);
  		free(weighted_C);
  	}
  	return model_;
  }
  
  void cross_validation(const problem *prob, const parameter *param, int nr_fold, double *target)
  {
  	int i;
  	int *fold_start;
  	int l = prob->l;
  	int *perm = Malloc(int,l);
  	if (nr_fold > l)
  	{
  		nr_fold = l;
  		fprintf(stderr,"WARNING: # folds > # data. Will use # folds = # data instead (i.e., leave-one-out cross validation)\n");
  	}
  	fold_start = Malloc(int,nr_fold+1);
  	for(i=0;i<l;i++) perm[i]=i;
  	for(i=0;i<l;i++)
  	{
  		int j = i+rand()%(l-i);
  		swap(perm[i],perm[j]);
  	}
  	for(i=0;i<=nr_fold;i++)
  		fold_start[i]=i*l/nr_fold;
  
  	for(i=0;i<nr_fold;i++)
  	{
  		int begin = fold_start[i];
  		int end = fold_start[i+1];
  		int j,k;
  		struct problem subprob;
  
  		subprob.bias = prob->bias;
  		subprob.n = prob->n;
  		subprob.l = l-(end-begin);
  		subprob.x = Malloc(struct feature_node*,subprob.l);
  		subprob.y = Malloc(double,subprob.l);
  
  		k=0;
  		for(j=0;j<begin;j++)
  		{
  			subprob.x[k] = prob->x[perm[j]];
  			subprob.y[k] = prob->y[perm[j]];
  			++k;
  		}
  		for(j=end;j<l;j++)
  		{
  			subprob.x[k] = prob->x[perm[j]];
  			subprob.y[k] = prob->y[perm[j]];
  			++k;
  		}
  		struct model *submodel = train(&subprob,param);
  		for(j=begin;j<end;j++)
  			target[perm[j]] = predict(submodel,prob->x[perm[j]]);
  		free_and_destroy_model(&submodel);
  		free(subprob.x);
  		free(subprob.y);
  	}
  	free(fold_start);
  	free(perm);
  }
  
  double predict_values(const struct model *model_, const struct feature_node *x, double *dec_values)
  {
  	int idx;
  	int n;
  	if(model_->bias>=0)
  		n=model_->nr_feature+1;
  	else
  		n=model_->nr_feature;
  	double *w=model_->w;
  	int nr_class=model_->nr_class;
  	int i;
  	int nr_w;
  	if(nr_class==2 && model_->param.solver_type != MCSVM_CS)
  		nr_w = 1;
  	else
  		nr_w = nr_class;
  
  	const feature_node *lx=x;
  	for(i=0;i<nr_w;i++)
  		dec_values[i] = 0;
  	for(; (idx=lx->index)!=-1; lx++)
  	{
  		// the dimension of testing data may exceed that of training
  		if(idx<=n)
  			for(i=0;i<nr_w;i++)
  				dec_values[i] += w[(idx-1)*nr_w+i]*lx->value;
  	}
  
  	if(nr_class==2)
  	{
  		if(model_->param.solver_type == L2R_L2LOSS_SVR ||
  		   model_->param.solver_type == L2R_L1LOSS_SVR_DUAL ||
  		   model_->param.solver_type == L2R_L2LOSS_SVR_DUAL)
  			return dec_values[0];
  		else
  			return (dec_values[0]>0)?model_->label[0]:model_->label[1];
  	}
  	else
  	{
  		int dec_max_idx = 0;
  		for(i=1;i<nr_class;i++)
  		{
  			if(dec_values[i] > dec_values[dec_max_idx])
  				dec_max_idx = i;
  		}
  		return model_->label[dec_max_idx];
  	}
  }
  
  double predict(const model *model_, const feature_node *x)
  {
  	double *dec_values = Malloc(double, model_->nr_class);
  	double label=predict_values(model_, x, dec_values);
  	free(dec_values);
  	return label;
  }
  
  double predict_probability(const struct model *model_, const struct feature_node *x, double* prob_estimates)
  {
  	if(check_probability_model(model_))
  	{
  		int i;
  		int nr_class=model_->nr_class;
  		int nr_w;
  		if(nr_class==2)
  			nr_w = 1;
  		else
  			nr_w = nr_class;
  
  		double label=predict_values(model_, x, prob_estimates);
  		for(i=0;i<nr_w;i++)
  			prob_estimates[i]=1/(1+exp(-prob_estimates[i]));
  
  		if(nr_class==2) // for binary classification
  			prob_estimates[1]=1.-prob_estimates[0];
  		else
  		{
  			double sum=0;
  			for(i=0; i<nr_class; i++)
  				sum+=prob_estimates[i];
  
  			for(i=0; i<nr_class; i++)
  				prob_estimates[i]=prob_estimates[i]/sum;
  		}
  
  		return label;
  	}
  	else
  		return 0;
  }
  
  static const char *solver_type_table[]=
  {
  	"L2R_LR", "L2R_L2LOSS_SVC_DUAL", "L2R_L2LOSS_SVC", "L2R_L1LOSS_SVC_DUAL", "MCSVM_CS",
  	"L1R_L2LOSS_SVC", "L1R_LR", "L2R_LR_DUAL",
  	"", "", "",
  	"L2R_L2LOSS_SVR", "L2R_L2LOSS_SVR_DUAL", "L2R_L1LOSS_SVR_DUAL", NULL
  };
  
  int save_model(const char *model_file_name, const struct model *model_)
  {
  	int i;
  	int nr_feature=model_->nr_feature;
  	int n;
  	const parameter& param = model_->param;
  
  	if(model_->bias>=0)
  		n=nr_feature+1;
  	else
  		n=nr_feature;
  	int w_size = n;
  	FILE *fp = fopen(model_file_name,"w");
  	if(fp==NULL) return -1;
  
  	char *old_locale = strdup(setlocale(LC_ALL, NULL));
  	setlocale(LC_ALL, "C");
  
  	int nr_w;
  	if(model_->nr_class==2 && model_->param.solver_type != MCSVM_CS)
  		nr_w=1;
  	else
  		nr_w=model_->nr_class;
  
  	fprintf(fp, "solver_type %s\n", solver_type_table[param.solver_type]);
  	fprintf(fp, "nr_class %d\n", model_->nr_class);
  
  	if(model_->label)
  	{
  		fprintf(fp, "label");
  		for(i=0; i<model_->nr_class; i++)
  			fprintf(fp, " %d", model_->label[i]);
  		fprintf(fp, "\n");
  	}
  
  	fprintf(fp, "nr_feature %d\n", nr_feature);
  
  	fprintf(fp, "bias %.16g\n", model_->bias);
  
  	fprintf(fp, "w\n");
  	for(i=0; i<w_size; i++)
  	{
  		int j;
  		for(j=0; j<nr_w; j++)
  			fprintf(fp, "%.16g ", model_->w[i*nr_w+j]);
  		fprintf(fp, "\n");
  	}
  
  	setlocale(LC_ALL, old_locale);
  	free(old_locale);
  
  	if (ferror(fp) != 0 || fclose(fp) != 0) return -1;
  	else return 0;
  }
  
  struct model *load_model(const char *model_file_name)
  {
  	FILE *fp = fopen(model_file_name,"r");
  	if(fp==NULL) return NULL;
  
  	int i;
  	int nr_feature;
  	int n;
  	int nr_class;
  	double bias;
  	model *model_ = Malloc(model,1);
  	parameter& param = model_->param;
  
  	model_->label = NULL;
  
  	char *old_locale = strdup(setlocale(LC_ALL, NULL));
  	setlocale(LC_ALL, "C");
  
  	char cmd[81];
  	while(1)
  	{
  		fscanf(fp,"%80s",cmd);
  		if(strcmp(cmd,"solver_type")==0)
  		{
  			fscanf(fp,"%80s",cmd);
  			int i;
  			for(i=0;solver_type_table[i];i++)
  			{
  				if(strcmp(solver_type_table[i],cmd)==0)
  				{
  					param.solver_type=i;
  					break;
  				}
  			}
  			if(solver_type_table[i] == NULL)
  			{
  				fprintf(stderr,"unknown solver type.\n");
  
  				setlocale(LC_ALL, old_locale);
  				free(model_->label);
  				free(model_);
  				free(old_locale);
  				return NULL;
  			}
  		}
  		else if(strcmp(cmd,"nr_class")==0)
  		{
  			fscanf(fp,"%d",&nr_class);
  			model_->nr_class=nr_class;
  		}
  		else if(strcmp(cmd,"nr_feature")==0)
  		{
  			fscanf(fp,"%d",&nr_feature);
  			model_->nr_feature=nr_feature;
  		}
  		else if(strcmp(cmd,"bias")==0)
  		{
  			fscanf(fp,"%lf",&bias);
  			model_->bias=bias;
  		}
  		else if(strcmp(cmd,"w")==0)
  		{
  			break;
  		}
  		else if(strcmp(cmd,"label")==0)
  		{
  			int nr_class = model_->nr_class;
  			model_->label = Malloc(int,nr_class);
  			for(int i=0;i<nr_class;i++)
  				fscanf(fp,"%d",&model_->label[i]);
  		}
  		else
  		{
  			fprintf(stderr,"unknown text in model file: [%s]\n",cmd);
  			setlocale(LC_ALL, old_locale);
  			free(model_->label);
  			free(model_);
  			free(old_locale);
  			return NULL;
  		}
  	}
  
  	nr_feature=model_->nr_feature;
  	if(model_->bias>=0)
  		n=nr_feature+1;
  	else
  		n=nr_feature;
  	int w_size = n;
  	int nr_w;
  	if(nr_class==2 && param.solver_type != MCSVM_CS)
  		nr_w = 1;
  	else
  		nr_w = nr_class;
  
  	model_->w=Malloc(double, w_size*nr_w);
  	for(i=0; i<w_size; i++)
  	{
  		int j;
  		for(j=0; j<nr_w; j++)
  			fscanf(fp, "%lf ", &model_->w[i*nr_w+j]);
  		fscanf(fp, "\n");
  	}
  
  	setlocale(LC_ALL, old_locale);
  	free(old_locale);
  
  	if (ferror(fp) != 0 || fclose(fp) != 0) return NULL;
  
  	return model_;
  }
  
  int get_nr_feature(const model *model_)
  {
  	return model_->nr_feature;
  }
  
  int get_nr_class(const model *model_)
  {
  	return model_->nr_class;
  }
  
  void get_labels(const model *model_, int* label)
  {
  	if (model_->label != NULL)
  		for(int i=0;i<model_->nr_class;i++)
  			label[i] = model_->label[i];
  }
  
  void free_model_content(struct model *model_ptr)
  {
  	if(model_ptr->w != NULL)
  		free(model_ptr->w);
  	if(model_ptr->label != NULL)
  		free(model_ptr->label);
  }
  
  void free_and_destroy_model(struct model **model_ptr_ptr)
  {
  	struct model *model_ptr = *model_ptr_ptr;
  	if(model_ptr != NULL)
  	{
  		free_model_content(model_ptr);
  		free(model_ptr);
  	}
  }
  
  void destroy_param(parameter* param)
  {
  	if(param->weight_label != NULL)
  		free(param->weight_label);
  	if(param->weight != NULL)
  		free(param->weight);
  }
  
  const char *check_parameter(const problem *prob, const parameter *param)
  {
  	if(param->eps <= 0)
  		return "eps <= 0";
  
  	if(param->C <= 0)
  		return "C <= 0";
  
  	if(param->p < 0)
  		return "p < 0";
  
  	if(param->solver_type != L2R_LR
  		&& param->solver_type != L2R_L2LOSS_SVC_DUAL
  		&& param->solver_type != L2R_L2LOSS_SVC
  		&& param->solver_type != L2R_L1LOSS_SVC_DUAL
  		&& param->solver_type != MCSVM_CS
  		&& param->solver_type != L1R_L2LOSS_SVC
  		&& param->solver_type != L1R_LR
  		&& param->solver_type != L2R_LR_DUAL
  		&& param->solver_type != L2R_L2LOSS_SVR
  		&& param->solver_type != L2R_L2LOSS_SVR_DUAL
  		&& param->solver_type != L2R_L1LOSS_SVR_DUAL)
  		return "unknown solver type";
  
  	return NULL;
  }
  
  int check_probability_model(const struct model *model_)
  {
  	return (model_->param.solver_type==L2R_LR ||
  			model_->param.solver_type==L2R_LR_DUAL ||
  			model_->param.solver_type==L1R_LR);
  }
  
  void set_print_string_function(void (*print_func)(const char*))
  {
  	if (print_func == NULL)
  		liblinear_print_string = &print_string_stdout;
  	else
  		liblinear_print_string = print_func;
  }
  
  