text: #include <math.h>
  #include <stdio.h>
  #include <stdlib.h>
  #include <ctype.h>
  #include <float.h>
  #include <string.h>
  #include <stdarg.h>
  #include <limits.h>
  #include <locale.h>
  #include <algorithm>
  #include "svm.h"
  int libsvm_version = LIBSVM_VERSION;
  typedef float Qfloat;
  typedef signed char schar;
  /*
  #ifndef min
  template <class T> static inline T std::min(T x,T y) { return (x<y)?x:y; }
  #endif
  #ifndef max
  template <class T> static inline T std::max(T x,T y) { return (x>y)?x:y; }
  #endif
  */
  //template <class T> static inline void std::swap(T& x, T& y) { T t=x; x=y; y=t; }
  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);
  }
  static inline double powi(double base, int times)
  {
  	double tmp = base, ret = 1.0;
  
  	for(int t=times; t>0; t/=2)
  	{
  		if(t%2==1) ret*=tmp;
  		tmp = tmp * tmp;
  	}
  	return ret;
  }
  #define INF HUGE_VAL
  #define TAU 1e-12
  #define Malloc(type,n) (type *)malloc((n)*sizeof(type))
  
  static void print_string_stdout(const char *s)
  {
  	fputs(s,stdout);
  	fflush(stdout);
  }
  static void (*svm_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);
  	(*svm_print_string)(buf);
  }
  #else
  static void info(const char *fmt,...) {}
  #endif
  
  //
  // Kernel Cache
  //
  // l is the number of total data items
  // size is the cache size limit in bytes
  //
  class Cache
  {
  public:
  	Cache(int l,long int size);
  	~Cache();
  
  	// request data [0,len)
  	// return some position p where [p,len) need to be filled
  	// (p >= len if nothing needs to be filled)
  	int get_data(const int index, Qfloat **data, int len);
  	void swap_index(int i, int j);	
  private:
  	int l;
  	long int size;
  	struct head_t
  	{
  		head_t *prev, *next;	// a circular list
  		Qfloat *data;
  		int len;		// data[0,len) is cached in this entry
  	};
  
  	head_t *head;
  	head_t lru_head;
  	void lru_delete(head_t *h);
  	void lru_insert(head_t *h);
  };
  
  Cache::Cache(int l_,long int size_):l(l_),size(size_)
  {
  	head = (head_t *)calloc(l,sizeof(head_t));	// initialized to 0
  	size /= sizeof(Qfloat);
  	size -= l * sizeof(head_t) / sizeof(Qfloat);
  	size = std::max(size, 2 * (long int) l);	// cache must be large enough for two columns
  	lru_head.next = lru_head.prev = &lru_head;
  }
  
  Cache::~Cache()
  {
  	for(head_t *h = lru_head.next; h != &lru_head; h=h->next)
  		free(h->data);
  	free(head);
  }
  
  void Cache::lru_delete(head_t *h)
  {
  	// delete from current location
  	h->prev->next = h->next;
  	h->next->prev = h->prev;
  }
  
  void Cache::lru_insert(head_t *h)
  {
  	// insert to last position
  	h->next = &lru_head;
  	h->prev = lru_head.prev;
  	h->prev->next = h;
  	h->next->prev = h;
  }
  
  int Cache::get_data(const int index, Qfloat **data, int len)
  {
  	head_t *h = &head[index];
  	if(h->len) lru_delete(h);
  	int more = len - h->len;
  
  	if(more > 0)
  	{
  		// free old space
  		while(size < more)
  		{
  			head_t *old = lru_head.next;
  			lru_delete(old);
  			free(old->data);
  			size += old->len;
  			old->data = 0;
  			old->len = 0;
  		}
  
  		// allocate new space
  		h->data = (Qfloat *)realloc(h->data,sizeof(Qfloat)*len);
  		size -= more;
  		std::swap(h->len,len);
  	}
  
  	lru_insert(h);
  	*data = h->data;
  	return len;
  }
  
  void Cache::swap_index(int i, int j)
  {
  	if(i==j) return;
  
  	if(head[i].len) lru_delete(&head[i]);
  	if(head[j].len) lru_delete(&head[j]);
  	std::swap(head[i].data,head[j].data);
  	std::swap(head[i].len,head[j].len);
  	if(head[i].len) lru_insert(&head[i]);
  	if(head[j].len) lru_insert(&head[j]);
  
  	if(i>j) std::swap(i,j);
  	for(head_t *h = lru_head.next; h!=&lru_head; h=h->next)
  	{
  		if(h->len > i)
  		{
  			if(h->len > j)
  				std::swap(h->data[i],h->data[j]);
  			else
  			{
  				// give up
  				lru_delete(h);
  				free(h->data);
  				size += h->len;
  				h->data = 0;
  				h->len = 0;
  			}
  		}
  	}
  }
  
  //
  // Kernel evaluation
  //
  // the static method k_function is for doing single kernel evaluation
  // the constructor of Kernel prepares to calculate the l*l kernel matrix
  // the member function get_Q is for getting one column from the Q Matrix
  //
  class QMatrix {
  public:
  	virtual Qfloat *get_Q(int column, int len) const = 0;
  	virtual double *get_QD() const = 0;
  	virtual void swap_index(int i, int j) const = 0;
  	virtual ~QMatrix() {}
  };
  
  class Kernel: public QMatrix {
  public:
  	Kernel(int l, svm_node * const * x, const svm_parameter& param);
  	virtual ~Kernel();
  
  	static double k_function(const svm_node *x, const svm_node *y,
  				 const svm_parameter& param);
  	virtual Qfloat *get_Q(int column, int len) const = 0;
  	virtual double *get_QD() const = 0;
  	virtual void swap_index(int i, int j) const	// no so const...
  	{
  		std::swap(x[i],x[j]);
  		if(x_square) std::swap(x_square[i],x_square[j]);
  	}
  protected:
  
  	double (Kernel::*kernel_function)(int i, int j) const;
  
  private:
  	const svm_node **x;
  	double *x_square;
  
  	// svm_parameter
  	const int kernel_type;
  	const int degree;
  	const double gamma;
  	const double coef0;
  
  	static double dot(const svm_node *px, const svm_node *py);
  	double kernel_linear(int i, int j) const
  	{
  		return dot(x[i],x[j]);
  	}
  	double kernel_poly(int i, int j) const
  	{
  		return powi(gamma*dot(x[i],x[j])+coef0,degree);
  	}
  	double kernel_rbf(int i, int j) const
  	{
  		return exp(-gamma*(x_square[i]+x_square[j]-2*dot(x[i],x[j])));
  	}
  	double kernel_sigmoid(int i, int j) const
  	{
  		return tanh(gamma*dot(x[i],x[j])+coef0);
  	}
  	double kernel_precomputed(int i, int j) const
  	{
  		return x[i][(int)(x[j][0].value)].value;
  	}
  };
  
  Kernel::Kernel(int l, svm_node * const * x_, const svm_parameter& param)
  :kernel_type(param.kernel_type), degree(param.degree),
   gamma(param.gamma), coef0(param.coef0)
  {
  	switch(kernel_type)
  	{
  		case LINEAR:
  			kernel_function = &Kernel::kernel_linear;
  			break;
  		case POLY:
  			kernel_function = &Kernel::kernel_poly;
  			break;
  		case RBF:
  			kernel_function = &Kernel::kernel_rbf;
  			break;
  		case SIGMOID:
  			kernel_function = &Kernel::kernel_sigmoid;
  			break;
  		case PRECOMPUTED:
  			kernel_function = &Kernel::kernel_precomputed;
  			break;
  	}
  
  	clone(x,x_,l);
  
  	if(kernel_type == RBF)
  	{
  		x_square = new double[l];
  		for(int i=0;i<l;i++)
  			x_square[i] = dot(x[i],x[i]);
  	}
  	else
  		x_square = 0;
  }
  
  Kernel::~Kernel()
  {
  	delete[] x;
  	delete[] x_square;
  }
  
  double Kernel::dot(const svm_node *px, const svm_node *py)
  {
  	double sum = 0;
  	while(px->index != -1 && py->index != -1)
  	{
  		if(px->index == py->index)
  		{
  			sum += px->value * py->value;
  			++px;
  			++py;
  		}
  		else
  		{
  			if(px->index > py->index)
  				++py;
  			else
  				++px;
  		}			
  	}
  	return sum;
  }
  
  double Kernel::k_function(const svm_node *x, const svm_node *y,
  			  const svm_parameter& param)
  {
  	switch(param.kernel_type)
  	{
  		case LINEAR:
  			return dot(x,y);
  		case POLY:
  			return powi(param.gamma*dot(x,y)+param.coef0,param.degree);
  		case RBF:
  		{
  			double sum = 0;
  			while(x->index != -1 && y->index !=-1)
  			{
  				if(x->index == y->index)
  				{
  					double d = x->value - y->value;
  					sum += d*d;
  					++x;
  					++y;
  				}
  				else
  				{
  					if(x->index > y->index)
  					{	
  						sum += y->value * y->value;
  						++y;
  					}
  					else
  					{
  						sum += x->value * x->value;
  						++x;
  					}
  				}
  			}
  
  			while(x->index != -1)
  			{
  				sum += x->value * x->value;
  				++x;
  			}
  
  			while(y->index != -1)
  			{
  				sum += y->value * y->value;
  				++y;
  			}
  			
  			return exp(-param.gamma*sum);
  		}
  		case SIGMOID:
  			return tanh(param.gamma*dot(x,y)+param.coef0);
  		case PRECOMPUTED:  //x: test (validation), y: SV
  			return x[(int)(y->value)].value;
  		default:
  			return 0;  // Unreachable 
  	}
  }
  
  // An SMO algorithm in Fan et al., JMLR 6(2005), p. 1889--1918
  // Solves:
  //
  //	min 0.5(\alpha^T Q \alpha) + p^T \alpha
  //
  //		y^T \alpha = \delta
  //		y_i = +1 or -1
  //		0 <= alpha_i <= Cp for y_i = 1
  //		0 <= alpha_i <= Cn for y_i = -1
  //
  // Given:
  //
  //	Q, p, y, Cp, Cn, and an initial feasible point \alpha
  //	l is the size of vectors and matrices
  //	eps is the stopping tolerance
  //
  // solution will be put in \alpha, objective value will be put in obj
  //
  class Solver {
  public:
  	Solver() {};
  	virtual ~Solver() {};
  
  	struct SolutionInfo {
  		double obj;
  		double rho;
  		double upper_bound_p;
  		double upper_bound_n;
  		double r;	// for Solver_NU
  	};
  
  	void Solve(int l, const QMatrix& Q, const double *p_, const schar *y_,
  		   double *alpha_, double Cp, double Cn, double eps,
  		   SolutionInfo* si, int shrinking);
  protected:
  	int active_size;
  	schar *y;
  	double *G;		// gradient of objective function
  	enum { LOWER_BOUND, UPPER_BOUND, FREE };
  	char *alpha_status;	// LOWER_BOUND, UPPER_BOUND, FREE
  	double *alpha;
  	const QMatrix *Q;
  	const double *QD;
  	double eps;
  	double Cp,Cn;
  	double *p;
  	int *active_set;
  	double *G_bar;		// gradient, if we treat free variables as 0
  	int l;
  	bool unshrink;	// XXX
  
  	double get_C(int i)
  	{
  		return (y[i] > 0)? Cp : Cn;
  	}
  	void update_alpha_status(int i)
  	{
  		if(alpha[i] >= get_C(i))
  			alpha_status[i] = UPPER_BOUND;
  		else if(alpha[i] <= 0)
  			alpha_status[i] = LOWER_BOUND;
  		else alpha_status[i] = FREE;
  	}
  	bool is_upper_bound(int i) { return alpha_status[i] == UPPER_BOUND; }
  	bool is_lower_bound(int i) { return alpha_status[i] == LOWER_BOUND; }
  	bool is_free(int i) { return alpha_status[i] == FREE; }
  	void swap_index(int i, int j);
  	void reconstruct_gradient();
  	virtual int select_working_set(int &i, int &j);
  	virtual double calculate_rho();
  	virtual void do_shrinking();
  private:
  	bool be_shrunk(int i, double Gmax1, double Gmax2);	
  };
  
  void Solver::swap_index(int i, int j)
  {
  	Q->swap_index(i,j);
  	std::swap(y[i],y[j]);
  	std::swap(G[i],G[j]);
  	std::swap(alpha_status[i],alpha_status[j]);
  	std::swap(alpha[i],alpha[j]);
  	std::swap(p[i],p[j]);
  	std::swap(active_set[i],active_set[j]);
  	std::swap(G_bar[i],G_bar[j]);
  }
  
  void Solver::reconstruct_gradient()
  {
  	// reconstruct inactive elements of G from G_bar and free variables
  
  	if(active_size == l) return;
  
  	int i,j;
  	int nr_free = 0;
  
  	for(j=active_size;j<l;j++)
  		G[j] = G_bar[j] + p[j];
  
  	for(j=0;j<active_size;j++)
  		if(is_free(j))
  			nr_free++;
  
  	if(2*nr_free < active_size)
  		info("\nWARNING: using -h 0 may be faster\n");
  
  	if (nr_free*l > 2*active_size*(l-active_size))
  	{
  		for(i=active_size;i<l;i++)
  		{
  			const Qfloat *Q_i = Q->get_Q(i,active_size);
  			for(j=0;j<active_size;j++)
  				if(is_free(j))
  					G[i] += alpha[j] * Q_i[j];
  		}
  	}
  	else
  	{
  		for(i=0;i<active_size;i++)
  			if(is_free(i))
  			{
  				const Qfloat *Q_i = Q->get_Q(i,l);
  				double alpha_i = alpha[i];
  				for(j=active_size;j<l;j++)
  					G[j] += alpha_i * Q_i[j];
  			}
  	}
  }
  
  void Solver::Solve(int l, const QMatrix& Q, const double *p_, const schar *y_,
  		   double *alpha_, double Cp, double Cn, double eps,
  		   SolutionInfo* si, int shrinking)
  {
  	this->l = l;
  	this->Q = &Q;
  	QD=Q.get_QD();
  	clone(p, p_,l);
  	clone(y, y_,l);
  	clone(alpha,alpha_,l);
  	this->Cp = Cp;
  	this->Cn = Cn;
  	this->eps = eps;
  	unshrink = false;
  
  	// initialize alpha_status
  	{
  		alpha_status = new char[l];
  		for(int i=0;i<l;i++)
  			update_alpha_status(i);
  	}
  
  	// initialize active set (for shrinking)
  	{
  		active_set = new int[l];
  		for(int i=0;i<l;i++)
  			active_set[i] = i;
  		active_size = l;
  	}
  
  	// initialize gradient
  	{
  		G = new double[l];
  		G_bar = new double[l];
  		int i;
  		for(i=0;i<l;i++)
  		{
  			G[i] = p[i];
  			G_bar[i] = 0;
  		}
  		for(i=0;i<l;i++)
  			if(!is_lower_bound(i))
  			{
  				const Qfloat *Q_i = Q.get_Q(i,l);
  				double alpha_i = alpha[i];
  				int j;
  				for(j=0;j<l;j++)
  					G[j] += alpha_i*Q_i[j];
  				if(is_upper_bound(i))
  					for(j=0;j<l;j++)
  						G_bar[j] += get_C(i) * Q_i[j];
  			}
  	}
  
  	// optimization step
  
  	int iter = 0;
  	int max_iter = std::max(10000000, l>INT_MAX/100 ? INT_MAX : 100*l);
  	int counter = std::min(l,1000)+1;
  	
  	while(iter < max_iter)
  	{
  		// show progress and do shrinking
  
  		if(--counter == 0)
  		{
  			counter = std::min(l,1000);
  			if(shrinking) do_shrinking();
  			info(".");
  		}
  
  		int i,j;
  		if(select_working_set(i,j)!=0)
  		{
  			// reconstruct the whole gradient
  			reconstruct_gradient();
  			// reset active set size and check
  			active_size = l;
  			info("*");
  			if(select_working_set(i,j)!=0)
  				break;
  			else
  				counter = 1;	// do shrinking next iteration
  		}
  		
  		++iter;
  
  		// update alpha[i] and alpha[j], handle bounds carefully
  		
  		const Qfloat *Q_i = Q.get_Q(i,active_size);
  		const Qfloat *Q_j = Q.get_Q(j,active_size);
  
  		double C_i = get_C(i);
  		double C_j = get_C(j);
  
  		double old_alpha_i = alpha[i];
  		double old_alpha_j = alpha[j];
  
  		if(y[i]!=y[j])
  		{
  			double quad_coef = QD[i]+QD[j]+2*Q_i[j];
  			if (quad_coef <= 0)
  				quad_coef = TAU;
  			double delta = (-G[i]-G[j])/quad_coef;
  			double diff = alpha[i] - alpha[j];
  			alpha[i] += delta;
  			alpha[j] += delta;
  			
  			if(diff > 0)
  			{
  				if(alpha[j] < 0)
  				{
  					alpha[j] = 0;
  					alpha[i] = diff;
  				}
  			}
  			else
  			{
  				if(alpha[i] < 0)
  				{
  					alpha[i] = 0;
  					alpha[j] = -diff;
  				}
  			}
  			if(diff > C_i - C_j)
  			{
  				if(alpha[i] > C_i)
  				{
  					alpha[i] = C_i;
  					alpha[j] = C_i - diff;
  				}
  			}
  			else
  			{
  				if(alpha[j] > C_j)
  				{
  					alpha[j] = C_j;
  					alpha[i] = C_j + diff;
  				}
  			}
  		}
  		else
  		{
  			double quad_coef = QD[i]+QD[j]-2*Q_i[j];
  			if (quad_coef <= 0)
  				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];
  			}
  		}
  	}
  
  	if(iter >= max_iter)
  	{
  		if(active_size < l)
  		{
  			// reconstruct the whole gradient to calculate objective value
  			reconstruct_gradient();
  			active_size = l;
  			info("*");
  		}
  		fprintf(stderr,"\nWARNING: reaching max number of iterations\n");
  	}
  
  	// 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] + p[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[] p;
  	delete[] y;
  	delete[] alpha;
  	delete[] alpha_status;
  	delete[] active_set;
  	delete[] G;
  	delete[] G_bar;
  }
  
  // return 1 if already optimal, return 0 otherwise
  int 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 = QD[i]+QD[j]-2.0*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 = QD[i]+QD[j]+2.0*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;
  }
  
  bool Solver::be_shrunk(int i, double Gmax1, double Gmax2)
  {
  	if(is_upper_bound(i))
  	{
  		if(y[i]==+1)
  			return(-G[i] > Gmax1);
  		else
  			return(-G[i] > Gmax2);
  	}
  	else if(is_lower_bound(i))
  	{
  		if(y[i]==+1)
  			return(G[i] > Gmax2);
  		else	
  			return(G[i] > Gmax1);
  	}
  	else
  		return(false);
  }
  
  void Solver::do_shrinking()
  {
  	int i;
  	double Gmax1 = -INF;		// max { -y_i * grad(f)_i | i in I_up(\alpha) }
  	double Gmax2 = -INF;		// max { y_i * grad(f)_i | i in I_low(\alpha) }
  
  	// find maximal violating pair first
  	for(i=0;i<active_size;i++)
  	{
  		if(y[i]==+1)	
  		{
  			if(!is_upper_bound(i))	
  			{
  				if(-G[i] >= Gmax1)
  					Gmax1 = -G[i];
  			}
  			if(!is_lower_bound(i))	
  			{
  				if(G[i] >= Gmax2)
  					Gmax2 = G[i];
  			}
  		}
  		else	
  		{
  			if(!is_upper_bound(i))	
  			{
  				if(-G[i] >= Gmax2)
  					Gmax2 = -G[i];
  			}
  			if(!is_lower_bound(i))	
  			{
  				if(G[i] >= Gmax1)
  					Gmax1 = G[i];
  			}
  		}
  	}
  
  	if(unshrink == false && Gmax1 + Gmax2 <= eps*10) 
  	{
  		unshrink = true;
  		reconstruct_gradient();
  		active_size = l;
  		info("*");
  	}
  
  	for(i=0;i<active_size;i++)
  		if (be_shrunk(i, Gmax1, Gmax2))
  		{
  			active_size--;
  			while (active_size > i)
  			{
  				if (!be_shrunk(active_size, Gmax1, Gmax2))
  				{
  					swap_index(i,active_size);
  					break;
  				}
  				active_size--;
  			}
  		}
  }
  
  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_upper_bound(i))
  		{
  			if(y[i]==-1)
  				ub = std::min(ub,yG);
  			else
  				lb = std::max(lb,yG);
  		}
  		else if(is_lower_bound(i))
  		{
  			if(y[i]==+1)
  				ub = std::min(ub,yG);
  			else
  				lb = std::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 *p, const schar *y,
  		   double *alpha, double Cp, double Cn, double eps,
  		   SolutionInfo* si, int shrinking)
  	{
  		this->si = si;
  		Solver::Solve(l,Q,p,y,alpha,Cp,Cn,eps,si,shrinking);
  	}
  private:
  	SolutionInfo *si;
  	int select_working_set(int &i, int &j);
  	double calculate_rho();
  	bool be_shrunk(int i, double Gmax1, double Gmax2, double Gmax3, double Gmax4);
  	void do_shrinking();
  };
  
  // return 1 if already optimal, return 0 otherwise
  int 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 = QD[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 = QD[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(std::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;
  }
  
  bool Solver_NU::be_shrunk(int i, double Gmax1, double Gmax2, double Gmax3, double Gmax4)
  {
  	if(is_upper_bound(i))
  	{
  		if(y[i]==+1)
  			return(-G[i] > Gmax1);
  		else	
  			return(-G[i] > Gmax4);
  	}
  	else if(is_lower_bound(i))
  	{
  		if(y[i]==+1)
  			return(G[i] > Gmax2);
  		else	
  			return(G[i] > Gmax3);
  	}
  	else
  		return(false);
  }
  
  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 i;
  	for(i=0;i<active_size;i++)
  	{
  		if(!is_upper_bound(i))
  		{
  			if(y[i]==+1)
  			{
  				if(-G[i] > Gmax1) Gmax1 = -G[i];
  			}
  			else	if(-G[i] > Gmax4) Gmax4 = -G[i];
  		}
  		if(!is_lower_bound(i))
  		{
  			if(y[i]==+1)
  			{	
  				if(G[i] > Gmax2) Gmax2 = G[i];
  			}
  			else	if(G[i] > Gmax3) Gmax3 = G[i];
  		}
  	}
  
  	if(unshrink == false && std::max(Gmax1+Gmax2,Gmax3+Gmax4) <= eps*10) 
  	{
  		unshrink = true;
  		reconstruct_gradient();
  		active_size = l;
  	}
  
  	for(i=0;i<active_size;i++)
  		if (be_shrunk(i, Gmax1, Gmax2, Gmax3, Gmax4))
  		{
  			active_size--;
  			while (active_size > i)
  			{
  				if (!be_shrunk(active_size, Gmax1, Gmax2, Gmax3, Gmax4))
  				{
  					swap_index(i,active_size);
  					break;
  				}
  				active_size--;
  			}
  		}
  }
  
  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_upper_bound(i))
  				lb1 = std::max(lb1,G[i]);
  			else if(is_lower_bound(i))
  				ub1 = std::min(ub1,G[i]);
  			else
  			{
  				++nr_free1;
  				sum_free1 += G[i];
  			}
  		}
  		else
  		{
  			if(is_upper_bound(i))
  				lb2 = std::max(lb2,G[i]);
  			else if(is_lower_bound(i))
  				ub2 = std::min(ub2,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,(long int)(param.cache_size*(1<<20)));
  		QD = new double[prob.l];
  		for(int i=0;i<prob.l;i++)
  			QD[i] = (this->*kernel_function)(i,i);
  	}
  	
  	Qfloat *get_Q(int i, int len) const
  	{
  		Qfloat *data;
  		int start, j;
  		if((start = cache->get_data(i,&data,len)) < len)
  		{
  			for(j=start;j<len;j++)
  				data[j] = (Qfloat)(y[i]*y[j]*(this->*kernel_function)(i,j));
  		}
  		return data;
  	}
  
  	double *get_QD() const
  	{
  		return QD;
  	}
  
  	void swap_index(int i, int j) const
  	{
  		cache->swap_index(i,j);
  		Kernel::swap_index(i,j);
  		std::swap(y[i],y[j]);
  		std::swap(QD[i],QD[j]);
  	}
  
  	~SVC_Q()
  	{
  		delete[] y;
  		delete cache;
  		delete[] QD;
  	}
  private:
  	schar *y;
  	Cache *cache;
  	double *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,(long int)(param.cache_size*(1<<20)));
  		QD = new double[prob.l];
  		for(int i=0;i<prob.l;i++)
  			QD[i] = (this->*kernel_function)(i,i);
  	}
  	
  	Qfloat *get_Q(int i, int len) const
  	{
  		Qfloat *data;
  		int start, j;
  		if((start = cache->get_data(i,&data,len)) < len)
  		{
  			for(j=start;j<len;j++)
  				data[j] = (Qfloat)(this->*kernel_function)(i,j);
  		}
  		return data;
  	}
  
  	double *get_QD() const
  	{
  		return QD;
  	}
  
  	void swap_index(int i, int j) const
  	{
  		cache->swap_index(i,j);
  		Kernel::swap_index(i,j);
  		std::swap(QD[i],QD[j]);
  	}
  
  	~ONE_CLASS_Q()
  	{
  		delete cache;
  		delete[] QD;
  	}
  private:
  	Cache *cache;
  	double *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,(long int)(param.cache_size*(1<<20)));
  		QD = new double[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] = (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
  	{
  		std::swap(sign[i],sign[j]);
  		std::swap(index[i],index[j]);
  		std::swap(QD[i],QD[j]);
  	}
  	
  	Qfloat *get_Q(int i, int len) const
  	{
  		Qfloat *data;
  		int j, real_i = index[i];
  		if(cache->get_data(real_i,&data,l) < l)
  		{
  			for(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(j=0;j<len;j++)
  			buf[j] = (Qfloat) si * (Qfloat) sign[j] * data[index[j]];
  		return buf;
  	}
  
  	double *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];
  	double *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] = std::min(1.0,sum_pos);
  			sum_pos -= alpha[i];
  		}
  		else
  		{
  			alpha[i] = std::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] = std::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;	
  };
  
  static 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;
  }
  
  // Platt's binary SVM Probablistic Output: an improvement from Lin et al.
  static 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-12;	// 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;
  			}
  			else
  				stepsize = stepsize / 2.0;
  		}
  
  		if (stepsize < min_step)
  		{
  			info("Line search fails in two-class probability estimates\n");
  			break;
  		}
  	}
  
  	if (iter>=max_iter)
  		info("Reaching maximal iterations in two-class probability estimates\n");
  	free(t);
  }
  
  static double sigmoid_predict(double decision_value, double A, double B)
  {
  	double fApB = decision_value*A+B;
  	// 1-p used later; avoid catastrophic cancellation
  	if (fApB >= 0)
  		return exp(-fApB)/(1.0+exp(-fApB));
  	else
  		return 1.0/(1+exp(fApB)) ;
  }
  
  // Method 2 from the multiclass_prob paper by Wu, Lin, and Weng
  static void multiclass_probability(int k, double **r, double *p)
  {
  	int t,j;
  	int iter = 0, max_iter=std::max(100,k);
  	double **Q=Malloc(double *,k);
  	double *Qp=Malloc(double,k);
  	double pQp, eps=0.005/k;
  	
  	for (t=0;t<k;t++)
  	{
  		p[t]=1.0/k;  // Valid if k = 1
  		Q[t]=Malloc(double,k);
  		Q[t][t]=0;
  		for (j=0;j<t;j++)
  		{
  			Q[t][t]+=r[j][t]*r[j][t];
  			Q[t][j]=Q[j][t];
  		}
  		for (j=t+1;j<k;j++)
  		{
  			Q[t][t]+=r[j][t]*r[j][t];
  			Q[t][j]=-r[j][t]*r[t][j];
  		}
  	}
  	for (iter=0;iter<max_iter;iter++)
  	{
  		// stopping condition, recalculate QP,pQP for numerical accuracy
  		pQp=0;
  		for (t=0;t<k;t++)
  		{
  			Qp[t]=0;
  			for (j=0;j<k;j++)
  				Qp[t]+=Q[t][j]*p[j];
  			pQp+=p[t]*Qp[t];
  		}
  		double max_error=0;
  		for (t=0;t<k;t++)
  		{
  			double error=fabs(Qp[t]-pQp);
  			if (error>max_error)
  				max_error=error;
  		}
  		if (max_error<eps) break;
  		
  		for (t=0;t<k;t++)
  		{
  			double diff=(-Qp[t]+pQp)/Q[t][t];
  			p[t]+=diff;
  			pQp=(pQp+diff*(diff*Q[t][t]+2*Qp[t]))/(1+diff)/(1+diff);
  			for (j=0;j<k;j++)
  			{
  				Qp[j]=(Qp[j]+diff*Q[t][j])/(1+diff);
  				p[j]/=(1+diff);
  			}
  		}
  	}
  	if (iter>=max_iter)
  		info("Exceeds max_iter in multiclass_prob\n");
  	for(t=0;t<k;t++) free(Q[t]);
  	free(Q);
  	free(Qp);
  }
  
  // Cross-validation decision values for probability estimates
  static void svm_binary_svc_probability(
  	const svm_problem *prob, const svm_parameter *param,
  	double Cp, double Cn, double& probA, double& probB)
  {
  	int i;
  	int nr_fold = 5;
  	int *perm = Malloc(int,prob->l);
  	double *dec_values = Malloc(double,prob->l);
  
  	// random shuffle
  	for(i=0;i<prob->l;i++) perm[i]=i;
  	for(i=0;i<prob->l;i++)
  	{
  		int j = i+rand()%(prob->l-i);
  		std::swap(perm[i],perm[j]);
  	}
  	for(i=0;i<nr_fold;i++)
  	{
  		int begin = i*prob->l/nr_fold;
  		int end = (i+1)*prob->l/nr_fold;
  		int j,k;
  		struct svm_problem subprob;
  
  		subprob.l = prob->l-(end-begin);
  		subprob.x = Malloc(struct svm_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<prob->l;j++)
  		{
  			subprob.x[k] = prob->x[perm[j]];
  			subprob.y[k] = prob->y[perm[j]];
  			++k;
  		}
  		int p_count=0,n_count=0;
  		for(j=0;j<k;j++)
  			if(subprob.y[j]>0)
  				p_count++;
  			else
  				n_count++;
  
  		if(p_count==0 && n_count==0)
  			for(j=begin;j<end;j++)
  				dec_values[perm[j]] = 0;
  		else if(p_count > 0 && n_count == 0)
  			for(j=begin;j<end;j++)
  				dec_values[perm[j]] = 1;
  		else if(p_count == 0 && n_count > 0)
  			for(j=begin;j<end;j++)
  				dec_values[perm[j]] = -1;
  		else
  		{
  			svm_parameter subparam = *param;
  			subparam.probability=0;
  			subparam.C=1.0;
  			subparam.nr_weight=2;
  			subparam.weight_label = Malloc(int,2);
  			subparam.weight = Malloc(double,2);
  			subparam.weight_label[0]=+1;
  			subparam.weight_label[1]=-1;
  			subparam.weight[0]=Cp;
  			subparam.weight[1]=Cn;
  			struct svm_model *submodel = svm_train(&subprob,&subparam);
  			for(j=begin;j<end;j++)
  			{
  				svm_predict_values(submodel,prob->x[perm[j]],&(dec_values[perm[j]])); 
  				// ensure +1 -1 order; reason not using CV subroutine
  				dec_values[perm[j]] *= submodel->label[0];
  			}		
  			svm_free_and_destroy_model(&submodel);
  			svm_destroy_param(&subparam);
  		}
  		free(subprob.x);
  		free(subprob.y);
  	}		
  	sigmoid_train(prob->l,dec_values,prob->y,probA,probB);
  	free(dec_values);
  	free(perm);
  }
  
  // Return parameter of a Laplace distribution 
  static double svm_svr_probability(
  	const svm_problem *prob, const svm_parameter *param)
  {
  	int i;
  	int nr_fold = 5;
  	double *ymv = Malloc(double,prob->l);
  	double mae = 0;
  
  	svm_parameter newparam = *param;
  	newparam.probability = 0;
  	svm_cross_validation(prob,&newparam,nr_fold,ymv);
  	for(i=0;i<prob->l;i++)
  	{
  		ymv[i]=prob->y[i]-ymv[i];
  		mae += fabs(ymv[i]);
  	}		
  	mae /= prob->l;
  	double std=sqrt(2*mae*mae);
  	int count=0;
  	mae=0;
  	for(i=0;i<prob->l;i++)
  		if (fabs(ymv[i]) > 5*std) 
  			count=count+1;
  		else 
  			mae+=fabs(ymv[i]);
  	mae /= (prob->l-count);
  	info("Prob. model for test data: target value = predicted value + z,\nz: Laplace distribution e^(-|z|/sigma)/(2sigma),sigma= %g\n",mae);
  	free(ymv);
  	return mae;
  }
  
  
  // 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 svm_group_classes(const svm_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)
  	{
  		std::swap(label[0],label[1]);
  		std::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);
  }
  
  //
  // Interface functions
  //
  svm_model *svm_train(const svm_problem *prob, const svm_parameter *param)
  {
  	svm_model *model = Malloc(svm_model,1);
  	model->param = *param;
  	model->free_sv = 0;	// XXX
  
  	if(param->svm_type == ONE_CLASS ||
  	   param->svm_type == EPSILON_SVR ||
  	   param->svm_type == NU_SVR)
  	{
  		// regression or one-class-svm
  		model->nr_class = 2;
  		model->label = NULL;
  		model->nSV = NULL;
  		model->probA = NULL; model->probB = NULL;
  		model->sv_coef = Malloc(double *,1);
  
  		if(param->probability && 
  		   (param->svm_type == EPSILON_SVR ||
  		    param->svm_type == NU_SVR))
  		{
  			model->probA = Malloc(double,1);
  			model->probA[0] = svm_svr_probability(prob,param);
  		}
  
  		decision_function f = svm_train_one(prob,param,0,0);
  		model->rho = Malloc(double,1);
  		model->rho[0] = f.rho;
  
  		int nSV = 0;
  		int i;
  		for(i=0;i<prob->l;i++)
  			if(fabs(f.alpha[i]) > 0) ++nSV;
  		model->l = nSV;
  		model->SV = Malloc(svm_node *,nSV);
  		model->sv_coef[0] = Malloc(double,nSV);
  		model->sv_indices = Malloc(int,nSV);
  		int j = 0;
  		for(i=0;i<prob->l;i++)
  			if(fabs(f.alpha[i]) > 0)
  			{
  				model->SV[j] = prob->x[i];
  				model->sv_coef[0][j] = f.alpha[i];
  				model->sv_indices[j] = i+1;
  				++j;
  			}		
  
  		free(f.alpha);
  	}
  	else
  	{
  		// classification
  		int l = prob->l;
  		int nr_class;
  		int *label = NULL;
  		int *start = NULL;
  		int *count = NULL;
  		int *perm = Malloc(int,l);
  
  		// group training data of the same class
  		svm_group_classes(prob,&nr_class,&label,&start,&count,perm);
  		if(nr_class == 1) 
  			info("WARNING: training data in only one class. See README for details.\n");
  		
  		svm_node **x = Malloc(svm_node *,l);
  		int i;
  		for(i=0;i<l;i++)
  			x[i] = prob->x[perm[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++)
  		{	
  			int j;
  			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];
  		}
  
  		// train k*(k-1)/2 models
  		
  		bool *nonzero = Malloc(bool,l);
  		for(i=0;i<l;i++)
  			nonzero[i] = false;
  		decision_function *f = Malloc(decision_function,nr_class*(nr_class-1)/2);
  
  		double *probA=NULL,*probB=NULL;
  		if (param->probability)
  		{
  			probA=Malloc(double,nr_class*(nr_class-1)/2);
  			probB=Malloc(double,nr_class*(nr_class-1)/2);
  		}
  
  		int p = 0;
  		for(i=0;i<nr_class;i++)
  			for(int j=i+1;j<nr_class;j++)
  			{
  				svm_problem sub_prob;
  				int si = start[i], sj = start[j];
  				int ci = count[i], cj = count[j];
  				sub_prob.l = ci+cj;
  				sub_prob.x = Malloc(svm_node *,sub_prob.l);
  				sub_prob.y = Malloc(double,sub_prob.l);
  				int k;
  				for(k=0;k<ci;k++)
  				{
  					sub_prob.x[k] = x[si+k];
  					sub_prob.y[k] = +1;
  				}
  				for(k=0;k<cj;k++)
  				{
  					sub_prob.x[ci+k] = x[sj+k];
  					sub_prob.y[ci+k] = -1;
  				}
  
  				if(param->probability)
  					svm_binary_svc_probability(&sub_prob,param,weighted_C[i],weighted_C[j],probA[p],probB[p]);
  
  				f[p] = svm_train_one(&sub_prob,param,weighted_C[i],weighted_C[j]);
  				for(k=0;k<ci;k++)
  					if(!nonzero[si+k] && fabs(f[p].alpha[k]) > 0)
  						nonzero[si+k] = true;
  				for(k=0;k<cj;k++)
  					if(!nonzero[sj+k] && fabs(f[p].alpha[ci+k]) > 0)
  						nonzero[sj+k] = true;
  				free(sub_prob.x);
  				free(sub_prob.y);
  				++p;
  			}
  
  		// build output
  
  		model->nr_class = nr_class;
  		
  		model->label = Malloc(int,nr_class);
  		for(i=0;i<nr_class;i++)
  			model->label[i] = label[i];
  		
  		model->rho = Malloc(double,nr_class*(nr_class-1)/2);
  		for(i=0;i<nr_class*(nr_class-1)/2;i++)
  			model->rho[i] = f[i].rho;
  
  		if(param->probability)
  		{
  			model->probA = Malloc(double,nr_class*(nr_class-1)/2);
  			model->probB = Malloc(double,nr_class*(nr_class-1)/2);
  			for(i=0;i<nr_class*(nr_class-1)/2;i++)
  			{
  				model->probA[i] = probA[i];
  				model->probB[i] = probB[i];
  			}
  		}
  		else
  		{
  			model->probA=NULL;
  			model->probB=NULL;
  		}
  
  		int total_sv = 0;
  		int *nz_count = Malloc(int,nr_class);
  		model->nSV = Malloc(int,nr_class);
  		for(i=0;i<nr_class;i++)
  		{
  			int nSV = 0;
  			for(int j=0;j<count[i];j++)
  				if(nonzero[start[i]+j])
  				{	
  					++nSV;
  					++total_sv;
  				}
  			model->nSV[i] = nSV;
  			nz_count[i] = nSV;
  		}
  		
  		info("Total nSV = %d\n",total_sv);
  
  		model->l = total_sv;
  		model->SV = Malloc(svm_node *,total_sv);
  		model->sv_indices = Malloc(int,total_sv);
  		p = 0;
  		for(i=0;i<l;i++)
  			if(nonzero[i])
  			{
  				model->SV[p] = x[i];
  				model->sv_indices[p++] = perm[i] + 1;
  			}
  
  		int *nz_start = Malloc(int,nr_class);
  		nz_start[0] = 0;
  		for(i=1;i<nr_class;i++)
  			nz_start[i] = nz_start[i-1]+nz_count[i-1];
  
  		model->sv_coef = Malloc(double *,nr_class-1);
  		for(i=0;i<nr_class-1;i++)
  			model->sv_coef[i] = Malloc(double,total_sv);
  
  		p = 0;
  		for(i=0;i<nr_class;i++)
  			for(int j=i+1;j<nr_class;j++)
  			{
  				// classifier (i,j): coefficients with
  				// i are in sv_coef[j-1][nz_start[i]...],
  				// j are in sv_coef[i][nz_start[j]...]
  
  				int si = start[i];
  				int sj = start[j];
  				int ci = count[i];
  				int cj = count[j];
  				
  				int q = nz_start[i];
  				int k;
  				for(k=0;k<ci;k++)
  					if(nonzero[si+k])
  						model->sv_coef[j-1][q++] = f[p].alpha[k];
  				q = nz_start[j];
  				for(k=0;k<cj;k++)
  					if(nonzero[sj+k])
  						model->sv_coef[i][q++] = f[p].alpha[ci+k];
  				++p;
  			}
  		
  		free(label);
  		free(probA);
  		free(probB);
  		free(count);
  		free(perm);
  		free(start);
  		free(x);
  		free(weighted_C);
  		free(nonzero);
  		for(i=0;i<nr_class*(nr_class-1)/2;i++)
  			free(f[i].alpha);
  		free(f);
  		free(nz_count);
  		free(nz_start);
  	}
  	return model;
  }
  
  // Stratified cross validation
  void svm_cross_validation(const svm_problem *prob, const svm_parameter *param, int nr_fold, double *target)
  {
  	int i;
  	int *fold_start;
  	int l = prob->l;
  	int *perm = Malloc(int,l);
  	int nr_class;
  	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);
  	// stratified cv may not give leave-one-out rate
  	// Each class to l folds -> some folds may have zero elements
  	if((param->svm_type == C_SVC ||
  	    param->svm_type == NU_SVC) && nr_fold < l)
  	{
  		int *start = NULL;
  		int *label = NULL;
  		int *count = NULL;
  		svm_group_classes(prob,&nr_class,&label,&start,&count,perm);
  
  		// random shuffle and then data grouped by fold using the array perm
  		int *fold_count = Malloc(int,nr_fold);
  		int c;
  		int *index = Malloc(int,l);
  		for(i=0;i<l;i++)
  			index[i]=perm[i];
  		for (c=0; c<nr_class; c++) 
  			for(i=0;i<count[c];i++)
  			{
  				int j = i+rand()%(count[c]-i);
  				std::swap(index[start[c]+j],index[start[c]+i]);
  			}
  		for(i=0;i<nr_fold;i++)
  		{
  			fold_count[i] = 0;
  			for (c=0; c<nr_class;c++)
  				fold_count[i]+=(i+1)*count[c]/nr_fold-i*count[c]/nr_fold;
  		}
  		fold_start[0]=0;
  		for (i=1;i<=nr_fold;i++)
  			fold_start[i] = fold_start[i-1]+fold_count[i-1];
  		for (c=0; c<nr_class;c++)
  			for(i=0;i<nr_fold;i++)
  			{
  				int begin = start[c]+i*count[c]/nr_fold;
  				int end = start[c]+(i+1)*count[c]/nr_fold;
  				for(int j=begin;j<end;j++)
  				{
  					perm[fold_start[i]] = index[j];
  					fold_start[i]++;
  				}
  			}
  		fold_start[0]=0;
  		for (i=1;i<=nr_fold;i++)
  			fold_start[i] = fold_start[i-1]+fold_count[i-1];
  		free(start);	
  		free(label);
  		free(count);	
  		free(index);
  		free(fold_count);
  	}
  	else
  	{
  		for(i=0;i<l;i++) perm[i]=i;
  		for(i=0;i<l;i++)
  		{
  			int j = i+rand()%(l-i);
  			std::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 svm_problem subprob;
  
  		subprob.l = l-(end-begin);
  		subprob.x = Malloc(struct svm_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 svm_model *submodel = svm_train(&subprob,param);
  		if(param->probability && 
  		   (param->svm_type == C_SVC || param->svm_type == NU_SVC))
  		{
  			double *prob_estimates=Malloc(double,svm_get_nr_class(submodel));
  			for(j=begin;j<end;j++)
  				target[perm[j]] = svm_predict_probability(submodel,prob->x[perm[j]],prob_estimates);
  			free(prob_estimates);			
  		}
  		else
  			for(j=begin;j<end;j++)
  				target[perm[j]] = svm_predict(submodel,prob->x[perm[j]]);
  		svm_free_and_destroy_model(&submodel);
  		free(subprob.x);
  		free(subprob.y);
  	}		
  	free(fold_start);
  	free(perm);	
  }
  
  
  int svm_get_svm_type(const svm_model *model)
  {
  	return model->param.svm_type;
  }
  
  int svm_get_nr_class(const svm_model *model)
  {
  	return model->nr_class;
  }
  
  void svm_get_labels(const svm_model *model, int* label)
  {
  	if (model->label != NULL)
  		for(int i=0;i<model->nr_class;i++)
  			label[i] = model->label[i];
  }
  
  void svm_get_sv_indices(const svm_model *model, int* indices)
  {
  	if (model->sv_indices != NULL)
  		for(int i=0;i<model->l;i++)
  			indices[i] = model->sv_indices[i];
  }
  
  int svm_get_nr_sv(const svm_model *model)
  {
  	return model->l;
  }
  
  double svm_get_svr_probability(const svm_model *model)
  {
  	if ((model->param.svm_type == EPSILON_SVR || model->param.svm_type == NU_SVR) &&
  	    model->probA!=NULL)
  		return model->probA[0];
  	else
  	{
  		fprintf(stderr,"Model doesn't contain information for SVR probability inference\n");
  		return 0;
  	}
  }
  
  double svm_predict_values(const svm_model *model, const svm_node *x, double* dec_values)
  {
  	int i;
  	if(model->param.svm_type == ONE_CLASS ||
  	   model->param.svm_type == EPSILON_SVR ||
  	   model->param.svm_type == NU_SVR)
  	{
  		double *sv_coef = model->sv_coef[0];
  		double sum = 0;
  		for(i=0;i<model->l;i++)
  			sum += sv_coef[i] * Kernel::k_function(x,model->SV[i],model->param);
  		sum -= model->rho[0];
  		*dec_values = sum;
  
  		if(model->param.svm_type == ONE_CLASS)
  			return (sum>0)?1:-1;
  		else
  			return sum;
  	}
  	else
  	{
  		int nr_class = model->nr_class;
  		int l = model->l;
  		
  		double *kvalue = Malloc(double,l);
  		for(i=0;i<l;i++)
  			kvalue[i] = Kernel::k_function(x,model->SV[i],model->param);
  
  		int *start = Malloc(int,nr_class);
  		start[0] = 0;
  		for(i=1;i<nr_class;i++)
  			start[i] = start[i-1]+model->nSV[i-1];
  
  		int *vote = Malloc(int,nr_class);
  		for(i=0;i<nr_class;i++)
  			vote[i] = 0;
  
  		int p=0;
  		for(i=0;i<nr_class;i++)
  			for(int j=i+1;j<nr_class;j++)
  			{
  				double sum = 0;
  				int si = start[i];
  				int sj = start[j];
  				int ci = model->nSV[i];
  				int cj = model->nSV[j];
  				
  				int k;
  				double *coef1 = model->sv_coef[j-1];
  				double *coef2 = model->sv_coef[i];
  				for(k=0;k<ci;k++)
  					sum += coef1[si+k] * kvalue[si+k];
  				for(k=0;k<cj;k++)
  					sum += coef2[sj+k] * kvalue[sj+k];
  				sum -= model->rho[p];
  				dec_values[p] = sum;
  
  				if(dec_values[p] > 0)
  					++vote[i];
  				else
  					++vote[j];
  				p++;
  			}
  
  		int vote_max_idx = 0;
  		for(i=1;i<nr_class;i++)
  			if(vote[i] > vote[vote_max_idx])
  				vote_max_idx = i;
  
  		free(kvalue);
  		free(start);
  		free(vote);
  		return model->label[vote_max_idx];
  	}
  }
  
  double svm_predict(const svm_model *model, const svm_node *x)
  {
  	int nr_class = model->nr_class;
  	double *dec_values;
  	if(model->param.svm_type == ONE_CLASS ||
  	   model->param.svm_type == EPSILON_SVR ||
  	   model->param.svm_type == NU_SVR)
  		dec_values = Malloc(double, 1);
  	else 
  		dec_values = Malloc(double, nr_class*(nr_class-1)/2);
  	double pred_result = svm_predict_values(model, x, dec_values);
  	free(dec_values);
  	return pred_result;
  }
  
  double svm_predict_probability(
  	const svm_model *model, const svm_node *x, double *prob_estimates)
  {
  	if ((model->param.svm_type == C_SVC || model->param.svm_type == NU_SVC) &&
  	    model->probA!=NULL && model->probB!=NULL)
  	{
  		int i;
  		int nr_class = model->nr_class;
  		double *dec_values = Malloc(double, nr_class*(nr_class-1)/2);
  		svm_predict_values(model, x, dec_values);
  
  		double min_prob=1e-7;
  		double **pairwise_prob=Malloc(double *,nr_class);
  		for(i=0;i<nr_class;i++)
  			pairwise_prob[i]=Malloc(double,nr_class);
  		int k=0;
  		for(i=0;i<nr_class;i++)
  			for(int j=i+1;j<nr_class;j++)
  			{
  				pairwise_prob[i][j]=std::min(std::max(sigmoid_predict(dec_values[k],model->probA[k],model->probB[k]),min_prob),1-min_prob);
  				pairwise_prob[j][i]=1-pairwise_prob[i][j];
  				k++;
  			}
  		multiclass_probability(nr_class,pairwise_prob,prob_estimates);
  
  		int prob_max_idx = 0;
  		for(i=1;i<nr_class;i++)
  			if(prob_estimates[i] > prob_estimates[prob_max_idx])
  				prob_max_idx = i;
  		for(i=0;i<nr_class;i++)
  			free(pairwise_prob[i]);
  		free(dec_values);
  		free(pairwise_prob);	     
  		return model->label[prob_max_idx];
  	}
  	else 
  		return svm_predict(model, x);
  }
  
  static const char *svm_type_table[] =
  {
  	"c_svc","nu_svc","one_class","epsilon_svr","nu_svr",NULL
  };
  
  static const char *kernel_type_table[]=
  {
  	"linear","polynomial","rbf","sigmoid","precomputed",NULL
  };
  
  int svm_save_model(const char *model_file_name, const svm_model *model)
  {
  	FILE *fp = fopen(model_file_name,"w");
  	if(fp==NULL) return -1;
  
  	char *old_locale = strdup(setlocale(LC_ALL, NULL));
  	setlocale(LC_ALL, "C");
  
  	const svm_parameter& param = model->param;
  
  	fprintf(fp,"svm_type %s\n", svm_type_table[param.svm_type]);
  	fprintf(fp,"kernel_type %s\n", kernel_type_table[param.kernel_type]);
  
  	if(param.kernel_type == POLY)
  		fprintf(fp,"degree %d\n", param.degree);
  
  	if(param.kernel_type == POLY || param.kernel_type == RBF || param.kernel_type == SIGMOID)
  		fprintf(fp,"gamma %g\n", param.gamma);
  
  	if(param.kernel_type == POLY || param.kernel_type == SIGMOID)
  		fprintf(fp,"coef0 %g\n", param.coef0);
  
  	int nr_class = model->nr_class;
  	int l = model->l;
  	fprintf(fp, "nr_class %d\n", nr_class);
  	fprintf(fp, "total_sv %d\n",l);
  	
  	{
  		fprintf(fp, "rho");
  		for(int i=0;i<nr_class*(nr_class-1)/2;i++)
  			fprintf(fp," %g",model->rho[i]);
  		fprintf(fp, "\n");
  	}
  	
  	if(model->label)
  	{
  		fprintf(fp, "label");
  		for(int i=0;i<nr_class;i++)
  			fprintf(fp," %d",model->label[i]);
  		fprintf(fp, "\n");
  	}
  
  	if(model->probA) // regression has probA only
  	{
  		fprintf(fp, "probA");
  		for(int i=0;i<nr_class*(nr_class-1)/2;i++)
  			fprintf(fp," %g",model->probA[i]);
  		fprintf(fp, "\n");
  	}
  	if(model->probB)
  	{
  		fprintf(fp, "probB");
  		for(int i=0;i<nr_class*(nr_class-1)/2;i++)
  			fprintf(fp," %g",model->probB[i]);
  		fprintf(fp, "\n");
  	}
  
  	if(model->nSV)
  	{
  		fprintf(fp, "nr_sv");
  		for(int i=0;i<nr_class;i++)
  			fprintf(fp," %d",model->nSV[i]);
  		fprintf(fp, "\n");
  	}
  
  	fprintf(fp, "SV\n");
  	const double * const *sv_coef = model->sv_coef;
  	const svm_node * const *SV = model->SV;
  
  	for(int i=0;i<l;i++)
  	{
  		for(int j=0;j<nr_class-1;j++)
  			fprintf(fp, "%.16g ",sv_coef[j][i]);
  
  		const svm_node *p = SV[i];
  
  		if(param.kernel_type == PRECOMPUTED)
  			fprintf(fp,"0:%d ",(int)(p->value));
  		else
  			while(p->index != -1)
  			{
  				fprintf(fp,"%d:%.8g ",p->index,p->value);
  				p++;
  			}
  		fprintf(fp, "\n");
  	}
  
  	setlocale(LC_ALL, old_locale);
  	free(old_locale);
  
  	if (ferror(fp) != 0 || fclose(fp) != 0) return -1;
  	else return 0;
  }
  
  static char *line = NULL;
  static int max_line_len;
  
  static char* readline(FILE *input)
  {
  	int len;
  
  	if(fgets(line,max_line_len,input) == NULL)
  		return NULL;
  
  	while(strrchr(line,'\n') == NULL)
  	{
  		max_line_len *= 2;
  		line = (char *) realloc(line,max_line_len);
  		len = (int) strlen(line);
  		if(fgets(line+len,max_line_len-len,input) == NULL)
  			break;
  	}
  	return line;
  }
  
  svm_model *svm_load_model(const char *model_file_name)
  {
  	FILE *fp = fopen(model_file_name,"rb");
  	if(fp==NULL) return NULL;
  
  	char *old_locale = strdup(setlocale(LC_ALL, NULL));
  	setlocale(LC_ALL, "C");
  
  	// read parameters
  
  	svm_model *model = Malloc(svm_model,1);
  	svm_parameter& param = model->param;
  	model->rho = NULL;
  	model->probA = NULL;
  	model->probB = NULL;
  	model->sv_indices = NULL;
  	model->label = NULL;
  	model->nSV = NULL;
  
  	char cmd[81];
  	while(1)
  	{
  		fscanf(fp,"%80s",cmd);
  
  		if(strcmp(cmd,"svm_type")==0)
  		{
  			fscanf(fp,"%80s",cmd);
  			int i;
  			for(i=0;svm_type_table[i];i++)
  			{
  				if(strcmp(svm_type_table[i],cmd)==0)
  				{
  					param.svm_type=i;
  					break;
  				}
  			}
  			if(svm_type_table[i] == NULL)
  			{
  				fprintf(stderr,"unknown svm type.\n");
  				
  				setlocale(LC_ALL, old_locale);
  				free(old_locale);
  				free(model->rho);
  				free(model->label);
  				free(model->nSV);
  				free(model);
  				return NULL;
  			}
  		}
  		else if(strcmp(cmd,"kernel_type")==0)
  		{		
  			fscanf(fp,"%80s",cmd);
  			int i;
  			for(i=0;kernel_type_table[i];i++)
  			{
  				if(strcmp(kernel_type_table[i],cmd)==0)
  				{
  					param.kernel_type=i;
  					break;
  				}
  			}
  			if(kernel_type_table[i] == NULL)
  			{
  				fprintf(stderr,"unknown kernel function.\n");
  				
  				setlocale(LC_ALL, old_locale);
  				free(old_locale);
  				free(model->rho);
  				free(model->label);
  				free(model->nSV);
  				free(model);
  				return NULL;
  			}
  		}
  		else if(strcmp(cmd,"degree")==0)
  			fscanf(fp,"%d",&param.degree);
  		else if(strcmp(cmd,"gamma")==0)
  			fscanf(fp,"%lf",&param.gamma);
  		else if(strcmp(cmd,"coef0")==0)
  			fscanf(fp,"%lf",&param.coef0);
  		else if(strcmp(cmd,"nr_class")==0)
  			fscanf(fp,"%d",&model->nr_class);
  		else if(strcmp(cmd,"total_sv")==0)
  			fscanf(fp,"%d",&model->l);
  		else if(strcmp(cmd,"rho")==0)
  		{
  			int n = model->nr_class * (model->nr_class-1)/2;
  			model->rho = Malloc(double,n);
  			for(int i=0;i<n;i++)
  				fscanf(fp,"%lf",&model->rho[i]);
  		}
  		else if(strcmp(cmd,"label")==0)
  		{
  			int n = model->nr_class;
  			model->label = Malloc(int,n);
  			for(int i=0;i<n;i++)
  				fscanf(fp,"%d",&model->label[i]);
  		}
  		else if(strcmp(cmd,"probA")==0)
  		{
  			int n = model->nr_class * (model->nr_class-1)/2;
  			model->probA = Malloc(double,n);
  			for(int i=0;i<n;i++)
  				fscanf(fp,"%lf",&model->probA[i]);
  		}
  		else if(strcmp(cmd,"probB")==0)
  		{
  			int n = model->nr_class * (model->nr_class-1)/2;
  			model->probB = Malloc(double,n);
  			for(int i=0;i<n;i++)
  				fscanf(fp,"%lf",&model->probB[i]);
  		}
  		else if(strcmp(cmd,"nr_sv")==0)
  		{
  			int n = model->nr_class;
  			model->nSV = Malloc(int,n);
  			for(int i=0;i<n;i++)
  				fscanf(fp,"%d",&model->nSV[i]);
  		}
  		else if(strcmp(cmd,"SV")==0)
  		{
  			while(1)
  			{
  				int c = getc(fp);
  				if(c==EOF || c=='\n') break;	
  			}
  			break;
  		}
  		else
  		{
  			fprintf(stderr,"unknown text in model file: [%s]\n",cmd);
  			
  			setlocale(LC_ALL, old_locale);
  			free(old_locale);
  			free(model->rho);
  			free(model->label);
  			free(model->nSV);
  			free(model);
  			return NULL;
  		}
  	}
  
  	// read sv_coef and SV
  
  	int elements = 0;
  	long pos = ftell(fp);
  
  	max_line_len = 1024;
  	line = Malloc(char,max_line_len);
  	char *p,*endptr,*idx,*val;
  
  	while(readline(fp)!=NULL)
  	{
  		p = strtok(line,":");
  		while(1)
  		{
  			p = strtok(NULL,":");
  			if(p == NULL)
  				break;
  			++elements;
  		}
  	}
  	elements += model->l;
  
  	fseek(fp,pos,SEEK_SET);
  
  	int m = model->nr_class - 1;
  	int l = model->l;
  	model->sv_coef = Malloc(double *,m);
  	int i;
  	for(i=0;i<m;i++)
  		model->sv_coef[i] = Malloc(double,l);
  	model->SV = Malloc(svm_node*,l);
  	svm_node *x_space = NULL;
  	if(l>0) x_space = Malloc(svm_node,elements);
  
  	int j=0;
  	for(i=0;i<l;i++)
  	{
  		readline(fp);
  		model->SV[i] = &x_space[j];
  
  		p = strtok(line, " \t");
  		model->sv_coef[0][i] = strtod(p,&endptr);
  		for(int k=1;k<m;k++)
  		{
  			p = strtok(NULL, " \t");
  			model->sv_coef[k][i] = strtod(p,&endptr);
  		}
  
  		while(1)
  		{
  			idx = strtok(NULL, ":");
  			val = strtok(NULL, " \t");
  
  			if(val == NULL)
  				break;
  			x_space[j].index = (int) strtol(idx,&endptr,10);
  			x_space[j].value = strtod(val,&endptr);
  
  			++j;
  		}
  		x_space[j++].index = -1;
  	}
  	free(line);
  
  	setlocale(LC_ALL, old_locale);
  	free(old_locale);
  
  	if (ferror(fp) != 0 || fclose(fp) != 0)
  		return NULL;
  
  	model->free_sv = 1;	// XXX
  	return model;
  }
  
  void svm_free_model_content(svm_model* model_ptr)
  {
  	if(model_ptr->free_sv && model_ptr->l > 0 && model_ptr->SV != NULL)
  		free((void *)(model_ptr->SV[0]));
  	if(model_ptr->sv_coef)
  	{
  		for(int i=0;i<model_ptr->nr_class-1;i++)
  			free(model_ptr->sv_coef[i]);
  	}
  
  	free(model_ptr->SV);
  	model_ptr->SV = NULL;
  
  	free(model_ptr->sv_coef);
  	model_ptr->sv_coef = NULL;
  
  	free(model_ptr->rho);
  	model_ptr->rho = NULL;
  
  	free(model_ptr->label);
  	model_ptr->label= NULL;
  
  	free(model_ptr->probA);
  	model_ptr->probA = NULL;
  
  	free(model_ptr->probB);
  	model_ptr->probB= NULL;
  
  	free(model_ptr->sv_indices);
  	model_ptr->sv_indices = NULL;
  
  	free(model_ptr->nSV);
  	model_ptr->nSV = NULL;
  }
  
  void svm_free_and_destroy_model(svm_model** model_ptr_ptr)
  {
  	if(model_ptr_ptr != NULL && *model_ptr_ptr != NULL)
  	{
  		svm_free_model_content(*model_ptr_ptr);
  		free(*model_ptr_ptr);
  		*model_ptr_ptr = NULL;
  	}
  }
  
  void svm_destroy_param(svm_parameter* param)
  {
  	free(param->weight_label);
  	free(param->weight);
  }
  
  const char *svm_check_parameter(const svm_problem *prob, const svm_parameter *param)
  {
  	// svm_type
  
  	int svm_type = param->svm_type;
  	if(svm_type != C_SVC &&
  	   svm_type != NU_SVC &&
  	   svm_type != ONE_CLASS &&
  	   svm_type != EPSILON_SVR &&
  	   svm_type != NU_SVR)
  		return "unknown svm type";
  	
  	// kernel_type, degree
  	
  	int kernel_type = param->kernel_type;
  	if(kernel_type != LINEAR &&
  	   kernel_type != POLY &&
  	   kernel_type != RBF &&
  	   kernel_type != SIGMOID &&
  	   kernel_type != PRECOMPUTED)
  		return "unknown kernel type";
  
  	if(param->gamma < 0)
  		return "gamma < 0";
  
  	if(param->degree < 0)
  		return "degree of polynomial kernel < 0";
  
  	// cache_size,eps,C,nu,p,shrinking
  
  	if(param->cache_size <= 0)
  		return "cache_size <= 0";
  
  	if(param->eps <= 0)
  		return "eps <= 0";
  
  	if(svm_type == C_SVC ||
  	   svm_type == EPSILON_SVR ||
  	   svm_type == NU_SVR)
  		if(param->C <= 0)
  			return "C <= 0";
  
  	if(svm_type == NU_SVC ||
  	   svm_type == ONE_CLASS ||
  	   svm_type == NU_SVR)
  		if(param->nu <= 0 || param->nu > 1)
  			return "nu <= 0 or nu > 1";
  
  	if(svm_type == EPSILON_SVR)
  		if(param->p < 0)
  			return "p < 0";
  
  	if(param->shrinking != 0 &&
  	   param->shrinking != 1)
  		return "shrinking != 0 and shrinking != 1";
  
  	if(param->probability != 0 &&
  	   param->probability != 1)
  		return "probability != 0 and probability != 1";
  
  	if(param->probability == 1 &&
  	   svm_type == ONE_CLASS)
  		return "one-class SVM probability output not supported yet";
  
  
  	// check whether nu-svc is feasible
  	
  	if(svm_type == NU_SVC)
  	{
  		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 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;
  				}
  			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;
  			}
  		}
  	
  		for(i=0;i<nr_class;i++)
  		{
  			int n1 = count[i];
  			for(int j=i+1;j<nr_class;j++)
  			{
  				int n2 = count[j];
  				if(param->nu*(n1+n2)/2 > std::min(n1,n2))
  				{
  					free(label);
  					free(count);
  					return "specified nu is infeasible";
  				}
  			}
  		}
  		free(label);
  		free(count);
  	}
  
  	return NULL;
  }
  
  int svm_check_probability_model(const svm_model *model)
  {
  	return ((model->param.svm_type == C_SVC || model->param.svm_type == NU_SVC) &&
  		model->probA!=NULL && model->probB!=NULL) ||
  		((model->param.svm_type == EPSILON_SVR || model->param.svm_type == NU_SVR) &&
  		 model->probA!=NULL);
  }
  
  void svm_set_print_string_function(void (*print_func)(const char *))
  {
  	if(print_func == NULL)
  		svm_print_string = &print_string_stdout;
  	else
  		svm_print_string = print_func;
  }
  