#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 "svm.h"

int libsvm_version = LIBSVM_VERSION;

typedef float Qfloat;
typedef signed char schar;

#ifndef min
template <class T>
static T min(const T x, const T y) { return (x < y) ? x : y; }
#endif

#ifndef max
template <class T>
static T max(const T x, const T y) { return (x > y) ? x : y; }
#endif

template <class T>
static void swap(T& x, T& y)
{
	T t = x;
	x   = y;
	y   = t;
}

template <class S, class T>
static void clone(T*& dst, S* src, const int n)
{
	dst = new T[n];
	memcpy((void*)dst, (void*)src, sizeof(T) * n);
}

static double powi(const double base, const 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 0
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(const int l_, const 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(const int l_, const 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          = max(size, 2 * long(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); }
	const 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 = nullptr;
			old->len  = 0;
		}

		// allocate new space
		h->data = (Qfloat*)realloc(h->data, sizeof(Qfloat) * len);
		size -= more;
		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]); }
	swap(head[i].data, head[j].data);
	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) { swap(i, j); }
	for (head_t* h = lru_head.next; h != &lru_head; h = h->next)
	{
		if (h->len > i)
		{
			if (h->len > j) { swap(h->data[i], h->data[j]); }
			else
			{
				// give up
				lru_delete(h);
				free(h->data);
				size += h->len;
				h->data = nullptr;
				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(const int l, svm_node* const* x_, const svm_parameter& param);
	~Kernel() override;

	static double k_function(const svm_node* x, const svm_node* y, const svm_parameter& param);
	Qfloat* get_Q(int column, int len) const override = 0;
	double* get_QD() const override = 0;
	void swap_index(const int i, const int j) const override
	// no so const...
	{
		swap(x[i], x[j]);
		if (x_square) { 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(const int i, const int j) const { return dot(x[i], x[j]); }
	double kernel_poly(const int i, const int j) const { return powi(gamma * dot(x[i], x[j]) + coef0, degree); }
	double kernel_rbf(const int i, const int j) const { return exp(-gamma * (x_square[i] + x_square[j] - 2 * dot(x[i], x[j]))); }
	double kernel_sigmoid(const int i, const int j) const { return tanh(gamma * dot(x[i], x[j]) + coef0); }
	double kernel_precomputed(const int i, const int j) const { return x[i][int(x[j][0].value)].value; }
};

Kernel::Kernel(const 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 = nullptr; }
}

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)
				{
					const 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:
	enum { LOWER_BOUND, UPPER_BOUND, FREE };

	int active_size = 0;
	schar* y        = nullptr;
	double* G       = nullptr;		// gradient of objective function

	char* alpha_status = nullptr;	// LOWER_BOUND, UPPER_BOUND, FREE
	double* alpha      = nullptr;
	const QMatrix* Q;
	const double* QD;
	double eps      = 0.0;
	double Cp       = 0.0, Cn = 0.0;
	double* p       = nullptr;
	int* active_set = nullptr;
	double* G_bar   = nullptr;		// gradient, if we treat free variables as 0
	int l           = 0.0;
	bool unshrink   = true;	// XXX

	double get_C(const int i) { return (y[i] > 0) ? Cp : Cn; }
	void update_alpha_status(const 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(const int i) { return alpha_status[i] == UPPER_BOUND; }
	bool is_lower_bound(const int i) { return alpha_status[i] == LOWER_BOUND; }
	bool is_free(const int i) { return alpha_status[i] == FREE; }
	void swap_index(const int i, const int j);
	void reconstruct_gradient();
	virtual int select_working_set(int& out_i, int& out_j);
	virtual double calculate_rho();
	virtual void do_shrinking();
private:
	bool be_shrunk(int i, double Gmax1, double Gmax2);
};

void Solver::swap_index(const int i, const int j)
{
	Q->swap_index(i, j);
	swap(y[i], y[j]);
	swap(G[i], G[j]);
	swap(alpha_status[i], alpha_status[j]);
	swap(alpha[i], alpha[j]);
	swap(p[i], p[j]);
	swap(active_set[i], active_set[j]);
	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);
				const double alpha_i = alpha[i];
				for (j = active_size; j < l; j++) { G[j] += alpha_i * Q_i[j]; }
			}
		}
	}
}

void Solver::Solve(const int l, const QMatrix& Q, const double* p_, const schar* y_, double* alpha_, const double Cp, const double Cn, const double eps,
				   SolutionInfo* si, const 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);
				const 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;
	const int max_iter = max(10000000, l > INT_MAX / 100 ? INT_MAX : 100 * l);
	int counter        = min(l, 1000) + 1;

	while (iter < max_iter)
	{
		// show progress and do shrinking

		if (--counter == 0)
		{
			counter = 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; } 	// do shrinking next iteration
			counter = 1;
		}

		++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);

		const double C_i = get_C(i);
		const double C_j = get_C(j);

		const double old_alpha_i = alpha[i];
		const 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; }
			const double delta = (-G[i] - G[j]) / quad_coef;
			const 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; }
			const double delta = (G[i] - G[j]) / quad_coef;
			const 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

		const double delta_alpha_i = alpha[i] - old_alpha_i;
		const 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

		{
			const bool ui = is_upper_bound(i);
			const 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;
		for (int 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;
				}
			}
		}
	}

	const int i       = Gmax_idx;
	const Qfloat* Q_i = nullptr;
	if (i != -1) { Q_i = Q->get_Q(i, active_size); } // NULL Q_i not accessed: Gmax=-INF if i=-1

	for (int j = 0; j < active_size; j++)
	{
		if (y[j] == +1)
		{
			if (!is_lower_bound(j))
			{
				const double grad_diff = Gmax + G[j];
				if (G[j] >= Gmax2) { Gmax2 = G[j]; }
				if (grad_diff > 0)
				{
					double obj_diff;
					const 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))
			{
				const double grad_diff = Gmax - G[j];
				if (-G[j] >= Gmax2) { Gmax2 = -G[j]; }
				if (grad_diff > 0)
				{
					double obj_diff;
					const 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 || Gmin_idx == -1) { return 1; }

	out_i = Gmax_idx;
	out_j = Gmin_idx;
	return 0;
}

bool Solver::be_shrunk(const int i, const double Gmax1, const double Gmax2)
{
	if (is_upper_bound(i))
	{
		if (y[i] == +1) { return (-G[i] > Gmax1); }
		return (-G[i] > Gmax2);
	}
	if (is_lower_bound(i))
	{
		if (y[i] == +1) { return (G[i] > Gmax2); }
		return (G[i] > Gmax1);
	}
	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++)
	{
		const double yG = y[i] * G[i];

		if (is_upper_bound(i))
		{
			if (y[i] == -1) { ub = min(ub, yG); }
			else { lb = max(lb, yG); }
		}
		else if (is_lower_bound(i))
		{
			if (y[i] == +1) { ub = min(ub, yG); }
			else { lb = max(lb, yG); }
		}
		else
		{
			++nr_free;
			sum_free += yG;
		}
	}

	if (nr_free > 0) { r = sum_free / nr_free; }
	else { r = (ub + lb) / 2; }

	return r;
}

//
// Solver for nu-svm classification and regression
//
// additional constraint: e^T \alpha = constant
//
class Solver_NU : public Solver
{
public:
	Solver_NU() {}
	void Solve(const int l, const QMatrix& Q, const double* p, const schar* y, double* alpha, const double Cp, const double Cn, const double eps,
			   SolutionInfo* si, const int shrinking)
	{
		this->si = si;
		Solver::Solve(l, Q, p, y, alpha, Cp, Cn, eps, si, shrinking);
	}
private:
	SolutionInfo* si = nullptr;
	int select_working_set(int& out_i, int& out_j) override;
	double calculate_rho() override;
	bool be_shrunk(const int i, const double Gmax1, const double Gmax2, const double Gmax3, const double Gmax4);
	void do_shrinking() override;
};

// 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;
				}
			}
		}
	}

	const int ip       = Gmaxp_idx;
	const int in       = Gmaxn_idx;
	const Qfloat* Q_ip = nullptr;
	const Qfloat* Q_in = nullptr;
	if (ip != -1) { Q_ip = Q->get_Q(ip, active_size); } // NULL Q_ip not accessed: Gmaxp=-INF if ip=-1
	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))
			{
				const double grad_diff = Gmaxp + G[j];
				if (G[j] >= Gmaxp2) { Gmaxp2 = G[j]; }
				if (grad_diff > 0)
				{
					double obj_diff;
					const 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))
			{
				const double grad_diff = Gmaxn - G[j];
				if (-G[j] >= Gmaxn2) { Gmaxn2 = -G[j]; }
				if (grad_diff > 0)
				{
					double obj_diff;
					const 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 (max(Gmaxp + Gmaxp2, Gmaxn + Gmaxn2) < eps || Gmin_idx == -1) { 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(const int i, const double Gmax1, const double Gmax2, const double Gmax3, const double Gmax4)
{
	if (is_upper_bound(i))
	{
		if (y[i] == +1) { return (-G[i] > Gmax1); }
		return (-G[i] > Gmax4);
	}
	if (is_lower_bound(i))
	{
		if (y[i] == +1) { return (G[i] > Gmax2); }
		return (G[i] > Gmax3);
	}
	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 && 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 = max(lb1, G[i]); }
			else if (is_lower_bound(i)) { ub1 = min(ub1, G[i]); }
			else
			{
				++nr_free1;
				sum_free1 += G[i];
			}
		}
		else
		{
			if (is_upper_bound(i)) { lb2 = max(lb2, G[i]); }
			else if (is_lower_bound(i)) { ub2 = 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(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(const int i, const int len) const override
	{
		Qfloat* data;
		int start;
		if ((start = cache->get_data(i, &data, len)) < len)
		{
			for (int j = start; j < len; j++) { data[j] = Qfloat(y[i] * y[j] * (this->*kernel_function)(i, j)); }
		}
		return data;
	}

	double* get_QD() const override { return QD; }

	void swap_index(const int i, const int j) const override
	{
		cache->swap_index(i, j);
		Kernel::swap_index(i, j);
		swap(y[i], y[j]);
		swap(QD[i], QD[j]);
	}

	~SVC_Q() override
	{
		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(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(const int i, const int len) const override
	{
		Qfloat* data;
		int start;
		if ((start = cache->get_data(i, &data, len)) < len) { for (int j = start; j < len; j++) { data[j] = Qfloat((this->*kernel_function)(i, j)); } }
		return data;
	}

	double* get_QD() const override { return QD; }

	void swap_index(const int i, const int j) const override
	{
		cache->swap_index(i, j);
		Kernel::swap_index(i, j);
		swap(QD[i], QD[j]);
	}

	~ONE_CLASS_Q() override
	{
		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(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(const int i, const int j) const override
	{
		swap(sign[i], sign[j]);
		swap(index[i], index[j]);
		swap(QD[i], QD[j]);
	}

	Qfloat* get_Q(const int i, const int len) const override
	{
		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;
		const 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 override { return QD; }

	~SVR_Q() override
	{
		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, const double Cp, const double Cn)
{
	const 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;
	const int l     = prob->l;
	const double nu = param->nu;

	schar* y = new schar[l];

	for (i = 0; i < l; i++)
	{
		if (prob->y[i] > 0) { y[i] = +1; }
		else { y[i] = -1; }
	}

	double sum_pos = nu * l / 2;
	double sum_neg = nu * l / 2;

	for (i = 0; i < l; i++)
	{
		if (y[i] == +1)
		{
			alpha[i] = min(1.0, sum_pos);
			sum_pos -= alpha[i];
		}
		else
		{
			alpha[i] = min(1.0, sum_neg);
			sum_neg -= alpha[i];
		}
	}

	double* zeros = new double[l];

	for (i = 0; i < l; i++) { zeros[i] = 0; }

	Solver_NU s;
	s.Solve(l, SVC_Q(*prob, *param, y), zeros, y, alpha, 1.0, 1.0, param->eps, si, param->shrinking);
	const 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)
{
	const int l   = prob->l;
	double* zeros = new double[l];
	schar* ones   = new schar[l];
	int i;

	const 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)
{
	const 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)
{
	const int l         = prob->l;
	const double C      = param->C;
	double* alpha2      = new double[2 * l];
	double* linear_term = new double[2 * l];
	schar* y            = new schar[2 * l];
	int i;

	double sum = C * param->nu * l / 2;
	for (i = 0; i < l; i++)
	{
		alpha2[i] = alpha2[i + l] = min(sum, C);
		sum -= alpha2[i];

		linear_term[i] = -prob->y[i];
		y[i]           = 1;

		linear_term[i + l] = prob->y[i];
		y[i + l]           = -1;
	}

	Solver_NU s;
	s.Solve(2 * l, SVR_Q(*prob, *param), linear_term, y, alpha2, C, C, param->eps, si, param->shrinking);

	info("epsilon = %f\n", -si->r);

	for (i = 0; i < l; i++) { alpha[i] = alpha2[i] - alpha2[i + l]; }

	delete[] alpha2;
	delete[] linear_term;
	delete[] y;
}

//
// decision_function
//
struct decision_function
{
	double* alpha;
	double rho;
};

static decision_function svm_train_one(const svm_problem* prob, const svm_parameter* param, const double Cp, const 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(const 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; }
	}

	const int max_iter    = 100;	// Maximal number of iterations
	const double min_step = 1e-10;	// Minimal step taken in line search
	const double sigma    = 1e-12;	// For numerically strict PD of Hessian
	const double eps      = 1e-5;
	const double hiTarget = (prior1 + 1.0) / (prior1 + 2.0);
	const double loTarget = 1 / (prior0 + 2.0);
	double* t             = Malloc(double, l);
	double fApB, p, q;
	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)
		double h11 = sigma; // numerically ensures strict PD
		double h22 = sigma;
		double h21 = 0.0;
		double g1  = 0.0;
		double 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));
			}
			const double d2 = p * q;
			h11 += dec_values[i] * dec_values[i] * d2;
			h22 += d2;
			h21 += dec_values[i] * d2;
			const double 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
		const double det = h11 * h22 - h21 * h21;
		const double dA  = -(h22 * g1 - h21 * g2) / det;
		const double dB  = -(-h21 * g1 + h11 * g2) / det;
		const double gd  = g1 * dA + g2 * dB;


		double stepsize = 1;		// Line Search
		while (stepsize >= min_step)
		{
			const double newA = A + stepsize * dA;
			const double newB = B + stepsize * dB;

			// New function value
			double 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;
			}
			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(const double decision_value, const double A, const double B)
{
	const double fApB = decision_value * A + B;
	// 1-p used later; avoid catastrophic cancellation
	if (fApB >= 0) { return exp(-fApB) / (1.0 + exp(-fApB)); }
	return 1.0 / (1 + exp(fApB));
}

// Method 2 from the multiclass_prob paper by Wu, Lin, and Weng
static void multiclass_probability(const int k, double** r, double* p)
{
	int t, j;
	int iter         = 0, max_iter = max(100, k);
	double** Q       = Malloc(double*, k);
	double* Qp       = Malloc(double, k);
	const double 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 < max_iter; iter++)
	{
		// stopping condition, recalculate QP,pQP for numerical accuracy
		double 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++)
		{
			const double error = fabs(Qp[t] - pQp);
			if (error > max_error) { max_error = error; }
		}
		if (max_error < eps) { break; }

		for (t = 0; t < k; t++)
		{
			const 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, const double Cp, const double Cn, double& probA, double& probB)
{
	int i;
	const 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++)
	{
		const int j = i + rand() % (prob->l - i);
		swap(perm[i], perm[j]);
	}
	for (i = 0; i < nr_fold; i++)
	{
		const int begin = i * prob->l / nr_fold;
		const int end   = (i + 1) * prob->l / nr_fold;
		int j;
		struct svm_problem subprob;

		subprob.l = prob->l - (end - begin);
		subprob.x = Malloc(struct svm_node*, subprob.l);
		subprob.y = Malloc(double, subprob.l);

		int 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;
	const 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;
	const 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)
{
	const 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++)
	{
		const 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);
}

//
// 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    = nullptr;
		model->nSV      = nullptr;
		model->probA    = nullptr;
		model->probB    = nullptr;
		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 = nullptr;
		int* start = nullptr;
		int* count = nullptr;
		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 = nullptr, *probB = nullptr;
		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 = nullptr;
			model->probB = nullptr;
		}

		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;
	const 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");
	}
	int* 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 = nullptr;
		int* label = nullptr;
		int* count = nullptr;
		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++)
			{
				const int j = i + rand() % (count[c] - i);
				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++)
			{
				const int begin = start[c] + i * count[c] / nr_fold;
				const 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++)
		{
			const 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++)
	{
		const int begin = fold_start[i];
		const int end   = fold_start[i + 1];
		int j;
		struct svm_problem subprob;

		subprob.l = l - (end - begin);
		subprob.x = Malloc(struct svm_node*, subprob.l);
		subprob.y = Malloc(double, subprob.l);

		int 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 != nullptr) { 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 != nullptr) { 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 != nullptr) { return model->probA[0]; }
	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; }
		return sum;
	}
	const int nr_class = model->nr_class;
	const 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;
			const int si = start[i];
			const int sj = start[j];
			const int ci = model->nSV[i];
			const 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)
{
	const 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); }
	const 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 != nullptr && model->probB != nullptr)
	{
		int i;
		const int nr_class = model->nr_class;
		double* dec_values = Malloc(double, nr_class * (nr_class - 1) / 2);
		svm_predict_values(model, x, dec_values);

		const 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] = min(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++;
			}
		}
		if (nr_class == 2)
		{
			prob_estimates[0] = pairwise_prob[0][1];
			prob_estimates[1] = pairwise_prob[1][0];
		}
		else { 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];
	}
	return svm_predict(model, x);
}

static const char* svm_type_table[]    = { "c_svc", "nu_svc", "one_class", "epsilon_svr", "nu_svr", nullptr };
static const char* kernel_type_table[] = { "linear", "polynomial", "rbf", "sigmoid", "precomputed", nullptr };

int svm_save_model(const char* model_file_name, const svm_model* model)
{
	FILE* fp = fopen(model_file_name, "w");
	if (fp == nullptr) { return -1; }

	char* old_locale = setlocale(LC_ALL, nullptr);
	if (old_locale) { old_locale = strdup(old_locale); }
	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 %.17g\n", param.gamma); }
	if (param.kernel_type == POLY || param.kernel_type == SIGMOID) { fprintf(fp, "coef0 %.17g\n", param.coef0); }

	const int nr_class = model->nr_class;
	const 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, " %.17g", 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, " %.17g", 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, " %.17g", 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, "%.17g ", 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; }
	return 0;
}

static char* line = nullptr;
static int max_line_len;

static char* readline(FILE* input)
{
	if (fgets(line, max_line_len, input) == nullptr) { return nullptr; }

	while (strrchr(line, '\n') == nullptr)
	{
		max_line_len *= 2;
		line          = (char*)realloc(line, max_line_len);
		const int len = int(strlen(line));
		if (fgets(line + len, max_line_len - len, input) == nullptr) { break; }
	}
	return line;
}

//
// FSCANF helps to handle fscanf failures.
// Its do-while block avoids the ambiguity when
// if (...)
//    FSCANF();
// is used
//
#define FSCANF(_stream, _format, _var) do{ if (fscanf(_stream, _format, _var) != 1) return false; }while(0)
bool read_model_header(FILE* fp, svm_model* model)
{
	svm_parameter& param = model->param;
	// parameters for training only won't be assigned, but arrays are assigned as NULL for safety
	param.nr_weight    = 0;
	param.weight_label = nullptr;
	param.weight       = nullptr;

	char cmd[81];
	while (true)
	{
		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] == nullptr)
			{
				fprintf(stderr, "unknown svm type.\n");
				return false;
			}
		}
		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] == nullptr)
			{
				fprintf(stderr, "unknown kernel function.\n");
				return false;
			}
		}
		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)
		{
			const 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)
		{
			const 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)
		{
			const 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)
		{
			const 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)
		{
			const 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 (true)
			{
				const int c = getc(fp);
				if (c == EOF || c == '\n') { break; }
			}
			break;
		}
		else
		{
			fprintf(stderr, "unknown text in model file: [%s]\n", cmd);
			return false;
		}
	}

	return true;
}

svm_model* svm_load_model(const char* model_file_name)
{
	FILE* fp = fopen(model_file_name, "rb");
	if (fp == nullptr) { return nullptr; }

	char* old_locale = setlocale(LC_ALL, nullptr);
	if (old_locale) { old_locale = strdup(old_locale); }
	setlocale(LC_ALL, "C");

	// read parameters

	svm_model* model  = Malloc(svm_model, 1);
	model->rho        = nullptr;
	model->probA      = nullptr;
	model->probB      = nullptr;
	model->sv_indices = nullptr;
	model->label      = nullptr;
	model->nSV        = nullptr;

	// read header
	if (!read_model_header(fp, model))
	{
		fprintf(stderr, "ERROR: fscanf failed to read model\n");
		setlocale(LC_ALL, old_locale);
		free(old_locale);
		free(model->rho);
		free(model->label);
		free(model->nSV);
		free(model);
		return nullptr;
	}

	// read sv_coef and SV

	int elements   = 0;
	const long pos = ftell(fp);

	max_line_len = 1024;
	line         = Malloc(char, max_line_len);
	char *p, *endptr;

	while (readline(fp) != nullptr)
	{
		p = strtok(line, ":");
		while (true)
		{
			p = strtok(nullptr, ":");
			if (p == nullptr) { break; }
			++elements;
		}
	}
	elements += model->l;

	fseek(fp, pos, SEEK_SET);

	const int m    = model->nr_class - 1;
	const 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 = nullptr;
	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(nullptr, " \t");
			model->sv_coef[k][i] = strtod(p, &endptr);
		}

		while (true)
		{
			char* idx = strtok(nullptr, ":");
			char* val = strtok(nullptr, " \t");

			if (val == nullptr) { 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 nullptr; }

	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 != nullptr) { 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 = nullptr;

	free(model_ptr->sv_coef);
	model_ptr->sv_coef = nullptr;

	free(model_ptr->rho);
	model_ptr->rho = nullptr;

	free(model_ptr->label);
	model_ptr->label = nullptr;

	free(model_ptr->probA);
	model_ptr->probA = nullptr;

	free(model_ptr->probB);
	model_ptr->probB = nullptr;

	free(model_ptr->sv_indices);
	model_ptr->sv_indices = nullptr;

	free(model_ptr->nSV);
	model_ptr->nSV = nullptr;
}

void svm_free_and_destroy_model(svm_model** model_ptr_ptr)
{
	if (model_ptr_ptr != nullptr && *model_ptr_ptr != nullptr)
	{
		svm_free_model_content(*model_ptr_ptr);
		free(*model_ptr_ptr);
		*model_ptr_ptr = nullptr;
	}
}

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

	const 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

	const 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 ((kernel_type == POLY || kernel_type == RBF || kernel_type == SIGMOID) && param->gamma < 0) { return "gamma < 0"; }
	if (kernel_type == POLY && 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)
	{
		const 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++)
		{
			const 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++)
		{
			const int n1 = count[i];
			for (int j = i + 1; j < nr_class; j++)
			{
				const int n2 = count[j];
				if (param->nu * (n1 + n2) / 2 > min(n1, n2))
				{
					free(label);
					free(count);
					return "specified nu is infeasible";
				}
			}
		}
		free(label);
		free(count);
	}

	return nullptr;
}

int svm_check_probability_model(const svm_model* model)
{
	return ((model->param.svm_type == C_SVC || model->param.svm_type == NU_SVC) && model->probA != nullptr && model->probB != nullptr)
		   || ((model->param.svm_type == EPSILON_SVR || model->param.svm_type == NU_SVR) && model->probA != nullptr);
}

void svm_set_print_string_function(void (*print_func)(const char*))
{
	if (print_func == nullptr) { svm_print_string = &print_string_stdout; }
	else { svm_print_string = print_func; }
}