CSvmTrainer.h

#ifndef SHARK_ALGORITHMS_CSVMTRAINER_H
#define SHARK_ALGORITHMS_CSVMTRAINER_H


#include <shark/Algorithms/Trainers/AbstractSvmTrainer.h>
#include <shark/Algorithms/Trainers/AbstractWeightedTrainer.h>
#include <shark/Algorithms/QP/BoxConstrainedProblems.h>
#include <shark/Algorithms/QP/SvmProblems.h>
#include <shark/Algorithms/QP/QpBoxLinear.h>
#include <shark/LinAlg/CachedMatrix.h>
#include <shark/LinAlg/GaussianKernelMatrix.h>
#include <shark/LinAlg/KernelMatrix.h>
#include <shark/LinAlg/PrecomputedMatrix.h>
#include <shark/LinAlg/RegularizedKernelMatrix.h>
#include <shark/Models/Kernels/GaussianRbfKernel.h>

//for MCSVMs!
#include <shark/Algorithms/QP/QpMcSimplexDecomp.h>
#include <shark/Algorithms/QP/QpMcBoxDecomp.h>
#include <shark/Algorithms/QP/QpMcLinear.h>
//~ #include <shark/Algorithms/Trainers/McSvm/McSvmMMRTrainer.h>
//~ #include <shark/Algorithms/Trainers/McSvm/McReinforcedSvmTrainer.h>

namespace shark {
	
	
enum class McSvm{
	WW,
	CS,
	LLW,
	ATM,
	ATS,
	ADM,
	OVA,
	MMR,
	ReinforcedSvm
};


///
/// \brief Training of C-SVMs for binary classification.
///
/// The C-SVM is the "standard" support vector machine for
/// binary (two-class) classification. Given are data tuples
/// \f$ (x_i, y_i) \f$ with x-component denoting input and
/// y-component denoting the label +1 or -1 (see the tutorial on
/// label conventions; the implementation uses values 0/1),
/// a kernel function k(x, x') and a regularization
/// constant C > 0. Let H denote the kernel induced
/// reproducing kernel Hilbert space of k, and let \f$ \phi \f$
/// denote the corresponding feature map.
/// Then the SVM classifier is the function
/// \f[
///     h(x) = \mathop{sign} (f(x))
/// \f]
/// \f[
///     f(x) = \langle w, \phi(x) \rangle + b
/// \f]
/// with coefficients w and b given by the (primal)
/// optimization problem
/// \f[
///     \min \frac{1}{2} \|w\|^2 + C \sum_i L(y_i, f(x_i)),
/// \f]
/// where \f$ L(y, f(x)) = \max\{0, 1 - y \cdot f(x)\} \f$
/// denotes the hinge loss.
///
/// For details refer to the paper:<br/>
/// <p>Support-Vector Networks. Corinna Cortes and Vladimir Vapnik,
/// Machine Learning, vol. 20 (1995), pp. 273-297.</p>
/// or simply to the Wikipedia article:<br/>
/// http://en.wikipedia.org/wiki/Support_vector_machine
/// \ingroup supervised_trainer
template <class InputType, class CacheType = float>
class CSvmTrainer : public AbstractSvmTrainer<
	InputType, unsigned int, 
	KernelClassifier<InputType>,
	AbstractWeightedTrainer<KernelClassifier<InputType> > 
>
{
private:
	typedef AbstractSvmTrainer<
		InputType, unsigned int, 
		KernelClassifier<InputType>,
		AbstractWeightedTrainer<KernelClassifier<InputType> > 
	> base_type;
public:

	/// \brief Convenience typedefs:
	/// this and many of the below typedefs build on the class template type CacheType.
	/// Simply changing that one template parameter CacheType thus allows to flexibly
	/// switch between using float or double as type for caching the kernel values.
	/// The default is float, offering sufficient accuracy in the vast majority
	/// of cases, at a memory cost of only four bytes. However, the template
	/// parameter makes it easy to use double instead, (e.g., in case high
	/// accuracy training is needed).
	typedef CacheType QpFloatType;

	typedef AbstractKernelFunction<InputType> KernelType;

	//! Constructor
	//! \param  kernel         kernel function to use for training and prediction
	//! \param  C              regularization parameter - always the 'true' value of C, even when unconstrained is set
	//! \param offset whether to train the svm with offset term
	//! \param  unconstrained  when a C-value is given via setParameter, should it be piped through the exp-function before using it in the solver?
	CSvmTrainer(KernelType* kernel, double C, bool offset, bool unconstrained = false)
	: base_type(kernel, C, offset, unconstrained), m_computeDerivative(false), m_McSvmType(McSvm::WW) //make  Vapnik happy!
	{ }
	
	//! Constructor
	//! \param  kernel         kernel function to use for training and prediction
	//! \param  negativeC   regularization parameter of the negative class (label 0)
	//! \param  positiveC    regularization parameter of the positive class (label 1)
	//! \param offset whether to train the svm with offset term
	//! \param  unconstrained  when a C-value is given via setParameter, should it be piped through the exp-function before using it in the solver?
	CSvmTrainer(KernelType* kernel, double negativeC, double positiveC, bool offset, bool unconstrained = false)
	: base_type(kernel,negativeC, positiveC, offset, unconstrained), m_computeDerivative(false), m_McSvmType(McSvm::WW) //make  Vapnik happy!
	{ }

	/// \brief From INameable: return the class name.
	std::string name() const
	{ return "CSvmTrainer"; }
	
	void setComputeBinaryDerivative(bool compute){
		m_computeDerivative = compute;
	}
	
	/// \brief sets the type of the multi-class svm used
	void setMcSvmType(McSvm type){
		m_McSvmType = type;
	}


	/// \brief Train the C-SVM.
	void train(KernelClassifier<InputType>& svm, LabeledData<InputType, unsigned int> const& dataset)
	{
		std::size_t classes = numberOfClasses(dataset);
		std::size_t ell = dataset.numberOfElements();
		if(classes == 2){
			// prepare model
			
			auto& f = svm.decisionFunction();
			if (f.basis() == dataset.inputs() && f.kernel() == base_type::m_kernel && f.alpha().size1() == ell && f.alpha().size2() == 1) {
				// warm start, keep the alphas (possibly clipped)
				if (this->m_trainOffset) f.offset() = RealVector(1);
			}
			else {
				f.setStructure(base_type::m_kernel, dataset.inputs(), this->m_trainOffset);
			}
			
			//dispatch to use the optimal implementation and solve the problem
			trainBinary(f,dataset);
			
			if (base_type::sparsify())
				f.sparsify();
			return;
		}
		
		//special case OVA 
		if(m_McSvmType == McSvm::OVA){
			trainOVA(svm,dataset);
			return;
		}
		
		//general multiclass case: find correct dual formulation
		bool sumToZero = false;
		bool simplex = false;
		QpSparseArray<QpFloatType> nu;
		QpSparseArray<QpFloatType> M;
		
		switch (m_McSvmType){
			case McSvm::WW:
				sumToZero = false;
				simplex = false;
				setupMcParametersWWCS(nu,M, classes);
			break;
			case McSvm::CS:
				sumToZero = false;
				simplex=true;
				setupMcParametersWWCS(nu,M, classes);
			break;
			case McSvm::LLW:
				sumToZero=true;
				simplex = false;
				setupMcParametersADMLLW(nu,M, classes);
			break;
			case McSvm::ATM:
				sumToZero=true;
				simplex=true;
				setupMcParametersATMATS(nu,M, classes);
			break;
			case McSvm::ATS:
				sumToZero=true;
				simplex = false;
				setupMcParametersATMATS(nu,M, classes);
			break;
			case McSvm::ADM:
				sumToZero=true;
				simplex=true;
				setupMcParametersADMLLW(nu,M, classes);
			break;
			case McSvm::ReinforcedSvm:
				sumToZero = false;
				simplex = false;
				setupMcParametersATMATS(nu,M, classes);
			break;
			case McSvm::MMR:
				sumToZero = true;
				simplex = true;
				setupMcParametersMMR(nu,M, classes);
			break;
			case McSvm::OVA: // handle OVA is switch statement to silence compiler warning
			break;
		}
		
		//setup linear part
		RealMatrix linear(ell,M.width(),1.0);
		if(m_McSvmType == McSvm::ReinforcedSvm){
			auto const& labels = dataset.labels();
			std::size_t i=0;
			for(unsigned int y: labels.elements()){
				linear(i, y) = classes - 1.0;   // self-margin target value of reinforced SVM loss
				i++;
			}
		}
		
		//solve dual
		RealMatrix alpha(ell,M.width(),0.0);
		RealVector bias(classes,0.0);
		if(simplex)
			solveMcSimplex(sumToZero,nu,M,linear,alpha,bias,dataset);
		else
			solveMcBox(sumToZero,nu,M,linear,alpha,bias,dataset);
		
		
		// write the solution into the model
		svm.decisionFunction().setStructure(this->m_kernel,dataset.inputs(),this->m_trainOffset,classes);
		
		for (std::size_t i=0; i<ell; i++)
		{
			unsigned int y = dataset.element(i).label;
			for (std::size_t c=0; c<classes; c++)
			{
				double sum = 0.0;
				std::size_t r = alpha.size2() * y;
				for (std::size_t p=0; p != alpha.size2(); p++, r++)
					sum += nu(r, c) * alpha(i, p);
				svm.decisionFunction().alpha(i,c) = sum;
			}
		}
		
		if (this->m_trainOffset) 
			svm.decisionFunction().offset() = bias;

		if (this->sparsify()) 
			svm.decisionFunction().sparsify();
	}

	/// \brief Train the C-SVM using weights.
	void train(KernelClassifier<InputType>& svm, WeightedLabeledData<InputType, unsigned int> const& dataset){
		SHARK_RUNTIME_CHECK(numberOfClasses(dataset) == 2, "CSVM with weights is only implemented for binary problems");
		// prepare model
		std::size_t n = dataset.numberOfElements();
		auto& f = svm.decisionFunction();
		if (f.basis() == dataset.inputs() && f.kernel() == base_type::m_kernel && f.alpha().size1() == n && f.alpha().size2() == 1) {
			// warm start, keep the alphas
			if (this->m_trainOffset) f.offset() = RealVector(1);
			else f.offset() = RealVector();
		}
		else {
			f.setStructure(base_type::m_kernel, dataset.inputs(), this->m_trainOffset);
		}

		//dispatch to use the optimal implementation and solve the problem
		trainBinary(f, dataset);

		if (base_type::sparsify()) f.sparsify();
	}
	
	RealVector const& get_db_dParams()const{
		return m_db_dParams;
	}

private:
	
	void solveMcSimplex(
		bool sumToZero, QpSparseArray<QpFloatType> const& nu,QpSparseArray<QpFloatType> const& M, RealMatrix const& linear,
		RealMatrix& alpha, RealVector& bias, 
		LabeledData<InputType, unsigned int> const& dataset
	){
		typedef KernelMatrix<InputType, QpFloatType> KernelMatrixType;
		typedef CachedMatrix< KernelMatrixType > CachedMatrixType;
		typedef PrecomputedMatrix< KernelMatrixType > PrecomputedMatrixType;
		
		KernelMatrixType km(*base_type::m_kernel, dataset.inputs());
		// solve the problem
		if (base_type::precomputeKernel())
		{
			PrecomputedMatrixType matrix(&km);
			QpMcSimplexDecomp< PrecomputedMatrixType> problem(matrix, M, dataset.labels(), linear, this->C());
			QpSolutionProperties& prop = base_type::m_solutionproperties;
			problem.setShrinking(base_type::m_shrinking);
			if(this->m_trainOffset){
				BiasSolverSimplex<PrecomputedMatrixType> biasSolver(&problem);
				biasSolver.solve(bias,base_type::m_stoppingcondition,nu,sumToZero, &prop);
			}
			else{
				QpSolver<QpMcSimplexDecomp< PrecomputedMatrixType> > solver(problem);
				solver.solve( base_type::m_stoppingcondition, &prop);
			}
			alpha = problem.solution();
		}
		else
		{
			CachedMatrixType matrix(&km, base_type::m_cacheSize);
			QpMcSimplexDecomp< CachedMatrixType> problem(matrix, M, dataset.labels(), linear, this->C());
			QpSolutionProperties& prop = base_type::m_solutionproperties;
			problem.setShrinking(base_type::m_shrinking);
			if(this->m_trainOffset){
				BiasSolverSimplex<CachedMatrixType> biasSolver(&problem);
				biasSolver.solve(bias,base_type::m_stoppingcondition,nu,sumToZero, &prop);
			}
			else{
				QpSolver<QpMcSimplexDecomp< CachedMatrixType> > solver(problem);
				solver.solve( base_type::m_stoppingcondition, &prop);
			}
			alpha = problem.solution();
		}
		base_type::m_accessCount = km.getAccessCount();
	}
	
	void solveMcBox(
		bool sumToZero, QpSparseArray<QpFloatType> const& nu,QpSparseArray<QpFloatType> const& M, RealMatrix const& linear,
		RealMatrix& alpha, RealVector& bias, 
		LabeledData<InputType, unsigned int> const& dataset
	){
		typedef KernelMatrix<InputType, QpFloatType> KernelMatrixType;
		typedef CachedMatrix< KernelMatrixType > CachedMatrixType;
		typedef PrecomputedMatrix< KernelMatrixType > PrecomputedMatrixType;
		
		KernelMatrixType km(*base_type::m_kernel, dataset.inputs());
		// solve the problem
		if (base_type::precomputeKernel())
		{
			PrecomputedMatrixType matrix(&km);
			QpMcBoxDecomp< PrecomputedMatrixType> problem(matrix, M, dataset.labels(), linear, this->C());
			QpSolutionProperties& prop = base_type::m_solutionproperties;
			problem.setShrinking(base_type::m_shrinking);
			if(this->m_trainOffset){
				BiasSolver<PrecomputedMatrixType> biasSolver(&problem);
				biasSolver.solve(bias,base_type::m_stoppingcondition,nu, sumToZero, &prop);
			}
			else{
				QpSolver<QpMcBoxDecomp< PrecomputedMatrixType> > solver(problem);
				solver.solve( base_type::m_stoppingcondition, &prop);
			}
			alpha = problem.solution();
		}
		else
		{
			CachedMatrixType matrix(&km, base_type::m_cacheSize);
			QpMcBoxDecomp< CachedMatrixType> problem(matrix, M, dataset.labels(), linear, this->C());
			QpSolutionProperties& prop = base_type::m_solutionproperties;
			problem.setShrinking(base_type::m_shrinking);
			if(this->m_trainOffset){
				BiasSolver<CachedMatrixType> biasSolver(&problem);
				biasSolver.solve(bias,base_type::m_stoppingcondition,nu, sumToZero, &prop);
			}
			else{
				QpSolver<QpMcBoxDecomp< CachedMatrixType> > solver(problem);
				solver.solve( base_type::m_stoppingcondition, &prop);
			}
			alpha = problem.solution();
		}
		base_type::m_accessCount = km.getAccessCount();
	}
	
	template<class Trainer>
	void trainMc(KernelClassifier<InputType>& svm, LabeledData<InputType, unsigned int> const& dataset){
		Trainer trainer(base_type::m_kernel,this->C(),this->m_trainOffset);
		trainer.stoppingCondition() = this->stoppingCondition();
		trainer.precomputeKernel() = this->precomputeKernel();
		trainer.sparsify() = this->sparsify();
		trainer.shrinking() = this->shrinking();
		trainer.s2do() = this->s2do();
		trainer.verbosity() = this->verbosity();
		trainer.setCacheSize(this->cacheSize());
		trainer.train(svm,dataset);
		this->solutionProperties() = trainer.solutionProperties();
		base_type::m_accessCount = trainer.accessCount();
	}
	
	void setupMcParametersWWCS(QpSparseArray<QpFloatType>& nu,QpSparseArray<QpFloatType>& M, std::size_t classes)const{
		nu.resize(classes * (classes-1), classes, 2*classes*(classes-1));
		for (unsigned int r=0, y=0; y<classes; y++)
		{
			for (unsigned int p=0, pp=0; p<classes-1; p++, pp++, r++)
			{
				if (pp == y) pp++;
				if (y < pp)
				{
					nu.add(r, y, 0.5);
					nu.add(r, pp, -0.5);
				}
				else
				{
					nu.add(r, pp, -0.5);
					nu.add(r, y, 0.5);
				}
			}
		}
		
		M.resize(classes * (classes-1) * classes, classes-1, 2 * classes * (classes-1) * (classes-1));
		for (unsigned int r=0, yv=0; yv<classes; yv++)
		{
			for (unsigned int pv=0, ppv=0; pv<classes-1; pv++, ppv++)
			{
				if (ppv == yv) ppv++;
				for (unsigned int yw=0; yw<classes; yw++, r++)
				{
					QpFloatType baseM = (yv == yw ? (QpFloatType)0.25 : (QpFloatType)0.0) - (ppv == yw ? (QpFloatType)0.25 : (QpFloatType)0.0);
					M.setDefaultValue(r, baseM);
					if (yv == yw)
					{
						M.add(r, ppv - (ppv >= yw ? 1 : 0), baseM + (QpFloatType)0.25);
					}
					else if (ppv == yw)
					{
						M.add(r, yv - (yv >= yw ? 1 : 0), baseM - (QpFloatType)0.25);
					}
					else
					{
						unsigned int pw = ppv - (ppv >= yw ? 1 : 0);
						unsigned int pw2 = yv - (yv >= yw ? 1 : 0);
						if (pw < pw2)
						{
							M.add(r, pw, baseM + (QpFloatType)0.25);
							M.add(r, pw2, baseM - (QpFloatType)0.25);
						}
						else
						{
							M.add(r, pw2, baseM - (QpFloatType)0.25);
							M.add(r, pw, baseM + (QpFloatType)0.25);
						}
					}
				}
			}
		}
	}
	
	void setupMcParametersATMATS(QpSparseArray<QpFloatType>& nu,QpSparseArray<QpFloatType>& M, std::size_t classes)const{
		nu.resize(classes*classes, classes, classes*classes);
		for (unsigned int r=0, y=0; y<classes; y++)
		{
			for (unsigned int p=0; p<classes; p++, r++)
			{
				nu.add(r, p, (QpFloatType)((p == y) ? 1.0 : -1.0));
			}
		}
		
		M.resize(classes * classes * classes, classes, 2 * classes * classes * classes);
		QpFloatType c_ne = (QpFloatType)(-1.0 / (double)classes);
		QpFloatType c_eq = (QpFloatType)1.0 + c_ne;
		for (unsigned int r=0, yv=0; yv<classes; yv++)
		{
			for (unsigned int pv=0; pv<classes; pv++)
			{
				QpFloatType sign = QpFloatType((yv == pv) ? -1 : 1);//cast to keep MSVC happy...
				for (unsigned int yw=0; yw<classes; yw++, r++)
				{
					M.setDefaultValue(r, sign * c_ne);
					if (yw == pv)
					{
						M.add(r, pv, -sign * c_eq);
					}
					else
					{
						M.add(r, pv, sign * c_eq);
						M.add(r, yw, -sign * c_ne);
					}
				}
			}
		}
	}
	
	void setupMcParametersADMLLW(QpSparseArray<QpFloatType>& nu,QpSparseArray<QpFloatType>& M, std::size_t classes)const{
		nu.resize(classes * (classes-1), classes, classes*(classes-1));
		for (unsigned int r=0, y=0; y<classes; y++)
		{
			for (unsigned int p=0, pp=0; p<classes-1; p++, pp++, r++)
			{
				if (pp == y) pp++;
				nu.add(r, pp, (QpFloatType)-1.0);
			}
		}
		
		M.resize(classes * (classes-1) * classes, classes-1, classes * (classes-1) * (classes-1));
		QpFloatType mood = (QpFloatType)(-1.0 / (double)classes);
		QpFloatType val = (QpFloatType)1.0 + mood;
		for (unsigned int r=0, yv=0; yv<classes; yv++)
		{
			for (unsigned int pv=0, ppv=0; pv<classes-1; pv++, ppv++)
			{
				if (ppv == yv) ppv++;
				for (unsigned int yw=0; yw<classes; yw++, r++)
				{
					M.setDefaultValue(r, mood);
					if (ppv != yw)
					{
						unsigned int pw = ppv - (ppv > yw ? 1 : 0);
						M.add(r, pw, val);
					}
				}
			}
		}
	}
	
	void setupMcParametersMMR(QpSparseArray<QpFloatType>& nu,QpSparseArray<QpFloatType>& M, std::size_t classes)const{
		nu.resize(classes, classes, classes);
		for (unsigned int y=0; y<classes; y++) 
			nu.add(y, y, 1.0);

		M.resize(classes * classes, 1, classes);
		QpFloatType mood = (QpFloatType)(-1.0 / (double)classes);
		QpFloatType val = (QpFloatType)1.0 + mood;
		for (unsigned int r=0, yv=0; yv<classes; yv++)
		{
			for (unsigned int yw=0; yw<classes; yw++, r++)
			{
				M.setDefaultValue(r, mood);
				if (yv == yw) M.add(r, 0, val);
			}
		}
	}
	
	void trainOVA(KernelClassifier<InputType>& svm, const LabeledData<InputType, unsigned int>& dataset){
		std::size_t classes = numberOfClasses(dataset);
		svm.decisionFunction().setStructure(this->m_kernel,dataset.inputs(),this->m_trainOffset,classes);
		
		base_type::m_solutionproperties.type = QpNone;
		base_type::m_solutionproperties.accuracy = 0.0;
		base_type::m_solutionproperties.iterations = 0;
		base_type::m_solutionproperties.value = 0.0;
		base_type::m_solutionproperties.seconds = 0.0;
		for (unsigned int c=0; c<classes; c++)
		{
			LabeledData<InputType, unsigned int> bindata = oneVersusRestProblem(dataset, c);
			KernelClassifier<InputType> binsvm;
// TODO: maybe build the Quadratic programs directly,
//       in order to profit from cached and
//       in particular from precomputed kernel
//       entries!
			CSvmTrainer<InputType, QpFloatType> bintrainer(base_type::m_kernel, this->C(),this->m_trainOffset);
			bintrainer.setCacheSize(this->cacheSize());
			bintrainer.sparsify() = false;
			bintrainer.stoppingCondition() = base_type::stoppingCondition();
			bintrainer.precomputeKernel() = base_type::precomputeKernel();		// sub-optimal!
			bintrainer.shrinking() = base_type::shrinking();
			bintrainer.s2do() = base_type::s2do();
			bintrainer.verbosity() = base_type::verbosity();
			bintrainer.train(binsvm, bindata);
			base_type::m_solutionproperties.iterations += bintrainer.solutionProperties().iterations;
			base_type::m_solutionproperties.seconds += bintrainer.solutionProperties().seconds;
			base_type::m_solutionproperties.accuracy = std::max(base_type::solutionProperties().accuracy, bintrainer.solutionProperties().accuracy);
			column(svm.decisionFunction().alpha(), c) = column(binsvm.decisionFunction().alpha(), 0);
			if (this->m_trainOffset)
				svm.decisionFunction().offset(c) = binsvm.decisionFunction().offset(0);
			base_type::m_accessCount += bintrainer.accessCount();
		}

		if (base_type::sparsify()) 
			svm.decisionFunction().sparsify();
	}
	
	//by default the normal unoptimized kernel matrix is used
	template<class T, class DatasetTypeT>
	void trainBinary(KernelExpansion<T>& svm, DatasetTypeT const& dataset){
		KernelMatrix<T, QpFloatType> km(*base_type::m_kernel, dataset.inputs());
		trainBinary(km,svm,dataset);
	}
	
	//in the case of a gaussian kernel and sparse vectors, we can use an optimized approach
	template<class T, class DatasetTypeT>
	void trainBinary(KernelExpansion<CompressedRealVector>& svm, DatasetTypeT const& dataset){
		//check whether a gaussian kernel is used
		typedef GaussianRbfKernel<CompressedRealVector> Gaussian;
		Gaussian const* kernel = dynamic_cast<Gaussian const*> (base_type::m_kernel);
		if(kernel != 0){//jep, use optimized kernel matrix
			GaussianKernelMatrix<CompressedRealVector,QpFloatType> km(kernel->gamma(),dataset.inputs());
			trainBinary(km,svm,dataset);
		}
		else{
			KernelMatrix<CompressedRealVector, QpFloatType> km(*base_type::m_kernel, dataset.inputs());
			trainBinary(km,svm,dataset);
		}
	}
	
	//create the problem for the unweighted datasets
	template<class Matrix, class T>
	void trainBinary(Matrix& km, KernelExpansion<T>& svm, LabeledData<T, unsigned int> const& dataset){
		if (QpConfig::precomputeKernel())
		{
			PrecomputedMatrix<Matrix> matrix(&km);
			CSVMProblem<PrecomputedMatrix<Matrix> > svmProblem(matrix,dataset.labels(),base_type::m_regularizers);
			optimize(svm,svmProblem,dataset);
		}
		else
		{
			CachedMatrix<Matrix> matrix(&km);
			CSVMProblem<CachedMatrix<Matrix> > svmProblem(matrix,dataset.labels(),base_type::m_regularizers);
			optimize(svm,svmProblem,dataset);
		}
		base_type::m_accessCount = km.getAccessCount();
	}
	
	// create the problem for the weighted datasets
	template<class Matrix, class T>
	void trainBinary(Matrix& km, KernelExpansion<T>& svm, WeightedLabeledData<T, unsigned int> const& dataset){
		if (QpConfig::precomputeKernel())
		{
			PrecomputedMatrix<Matrix> matrix(&km);
			GeneralQuadraticProblem<PrecomputedMatrix<Matrix> > svmProblem(
				matrix,dataset.labels(),dataset.weights(),base_type::m_regularizers
			);
			optimize(svm,svmProblem,dataset.data());
		}
		else
		{
			CachedMatrix<Matrix> matrix(&km);
			GeneralQuadraticProblem<CachedMatrix<Matrix> > svmProblem(
				matrix,dataset.labels(),dataset.weights(),base_type::m_regularizers
			);
			optimize(svm,svmProblem,dataset.data());
		}
		base_type::m_accessCount = km.getAccessCount();
	}
	
	template<class SVMProblemType>
	void optimize(KernelExpansion<InputType>& svm, SVMProblemType& svmProblem, LabeledData<InputType, unsigned int> const& dataset){
		if (this->m_trainOffset)
		{
			typedef SvmShrinkingProblem<SVMProblemType> ProblemType;
			ProblemType problem(svmProblem,base_type::m_shrinking);
			QpSolver< ProblemType > solver(problem);
			// truncate the existing solution to the bounds
			RealVector const& reg = this->regularizationParameters();
			double C_minus = reg(0);
			double C_plus = (reg.size() == 1) ? reg(0) : reg(1);
			std::size_t i=0;
			for (auto label : dataset.labels().elements()) {
				double a = svm.alpha()(i, 0);
				if (label == 0) a = std::max(std::min(a, 0.0), -C_minus);
				else            a = std::min(std::max(a, 0.0), C_plus);
				svm.alpha()(i, 0) = a;
				i++;
			}
			problem.setInitialSolution(blas::column(svm.alpha(), 0));
			solver.solve(base_type::stoppingCondition(), &base_type::solutionProperties());
			column(svm.alpha(),0)= problem.getUnpermutedAlpha();
			svm.offset(0) = computeBias(problem,dataset);
		}
		else
		{
			typedef BoxConstrainedShrinkingProblem<SVMProblemType> ProblemType;
			ProblemType problem(svmProblem,base_type::m_shrinking);
			QpSolver< ProblemType> solver(problem);
			// truncate the existing solution to the bounds
			RealVector const& reg = this->regularizationParameters();
			double C_minus = reg(0);
			double C_plus = (reg.size() == 1) ? reg(0) : reg(1);
			std::size_t i=0;
			for (auto label : dataset.labels().elements()) {
				double a = svm.alpha()(i, 0);
				if (label == 0) a = std::max(std::min(a, 0.0), -C_minus);
				else            a = std::min(std::max(a, 0.0), C_plus);
				svm.alpha()(i, 0) = a;
				i++;
			}
			problem.setInitialSolution(blas::column(svm.alpha(), 0));
			solver.solve(base_type::stoppingCondition(), &base_type::solutionProperties());
			column(svm.alpha(),0) = problem.getUnpermutedAlpha();
		}
	}
	
	RealVector m_db_dParams; ///< in the rare case that there are only bounded SVs and no free SVs, this will hold the derivative of b w.r.t. the hyperparameters. Derivative w.r.t. C is last.

	bool m_computeDerivative;
	McSvm m_McSvmType;

	template<class Problem>
	double computeBias(Problem const& problem, LabeledData<InputType, unsigned int> const& dataset){
		std::size_t nkp = base_type::m_kernel->numberOfParameters();
		m_db_dParams.resize(nkp+1);
		m_db_dParams.clear();

		std::size_t ell = problem.dimensions();
		if (ell == 0) return 0.0;

		// compute the offset from the KKT conditions
		double lowerBound = -1e100;
		double upperBound = 1e100;
		double sum = 0.0;
		std::size_t freeVars = 0;
		std::size_t lower_i = 0;
		std::size_t upper_i = 0;
		for (std::size_t i=0; i<ell; i++)
		{
			double value = problem.gradient(i);
			if (problem.alpha(i) == problem.boxMin(i))
			{
				if (value > lowerBound) { //in case of no free SVs, we are looking for the largest gradient of all alphas at the lower bound
					lowerBound = value;
					lower_i = i;
				}
			}
			else if (problem.alpha(i) == problem.boxMax(i))
			{
				if (value < upperBound) { //in case of no free SVs, we are looking for the smallest gradient of all alphas at the upper bound
					upperBound = value;
					upper_i = i;
				}
			}
			else
			{
				sum += value;
				freeVars++;
			}
		}
		if (freeVars > 0)
			return sum / freeVars;		//stabilized (averaged) exact value

		if(!m_computeDerivative)
			return 0.5 * (lowerBound + upperBound);	//best estimate

		lower_i = problem.permutation(lower_i);
		upper_i = problem.permutation(upper_i);

		SHARK_RUNTIME_CHECK(base_type::m_regularizers.size() == 1, "derivative only implemented for SVM with one C" );

		// We next compute the derivative of lowerBound and upperBound wrt C, in order to then get that of b wrt C.
		// The equation at the foundation of this simply is g_i = y_i - \sum_j \alpha_j K_{ij} .
		double dlower_dC = 0.0;
		double dupper_dC = 0.0;
		// At the same time, we also compute the derivative of lowerBound and upperBound wrt the kernel parameters.
		// The equation at the foundation of this simply is g_i = y_i - \sum_j \alpha_j K_{ij} .
		RealVector dupper_dkernel( nkp,0 );
		RealVector dlower_dkernel( nkp,0 );
		//state for eval and evalDerivative of the kernel
		boost::shared_ptr<State> kernelState = base_type::m_kernel->createState();
		RealVector der(nkp ); //derivative storage helper
		//todo: O.K.: here kernel single input derivative would be usefull
		//also it can be usefull to use here real batch processing and use batches of size 1 for lower /upper
		//and instead of singleInput whole batches.
		//what we do is, that we use the batched input versions with batches of size one.
		typename Batch<InputType>::type singleInput = Batch<InputType>::createBatch( dataset.element(0).input, 1 );
		typename Batch<InputType>::type lowerInput = Batch<InputType>::createBatch( dataset.element(lower_i).input, 1 );
		typename Batch<InputType>::type upperInput = Batch<InputType>::createBatch( dataset.element(upper_i).input, 1 );
		getBatchElement( lowerInput, 0 ) = dataset.element(lower_i).input; //copy the current input into the batch
		getBatchElement( upperInput, 0 ) = dataset.element(upper_i).input; //copy the current input into the batch
		RealMatrix one(1,1,1); //weight of input
		RealMatrix result(1,1); //stores the result of the call

		for (std::size_t i=0; i<ell; i++) {
			double cur_alpha = problem.alpha(problem.permutation(i));
			if ( cur_alpha != 0 ) {
				int cur_label = ( cur_alpha>0.0 ? 1 : -1 );
				getBatchElement( singleInput, 0 ) = dataset.element(i).input; //copy the current input into the batch
				// treat contributions of largest gradient at lower bound
				base_type::m_kernel->eval( lowerInput, singleInput, result, *kernelState );
				dlower_dC += cur_label * result(0,0);
				base_type::m_kernel->weightedParameterDerivative( lowerInput, singleInput,one, *kernelState, der );
				for ( std::size_t k=0; k<nkp; k++ ) {
					dlower_dkernel(k) += cur_label * der(k);
				}
				// treat contributions of smallest gradient at upper bound
				base_type::m_kernel->eval( upperInput, singleInput,result, *kernelState );
				dupper_dC += cur_label * result(0,0);
				base_type::m_kernel->weightedParameterDerivative( upperInput, singleInput, one, *kernelState, der );
				for ( std::size_t k=0; k<nkp; k++ ) {
					dupper_dkernel(k) += cur_label * der(k);
				}
			}
		}
		// assign final values to derivative of b wrt hyperparameters
		m_db_dParams( nkp ) = -0.5 * ( dlower_dC + dupper_dC );
		for ( std::size_t k=0; k<nkp; k++ ) {
			m_db_dParams(k) = -0.5 * this->C() * ( dlower_dkernel(k) + dupper_dkernel(k) );
		}
		if ( base_type::m_unconstrained ) {
			m_db_dParams( nkp ) *= this->C();
		}
		
		return 0.5 * (lowerBound + upperBound);	//best estimate
	}
};


template <class InputType>
class LinearCSvmTrainer : public AbstractLinearSvmTrainer<InputType>
{
public:
	typedef AbstractLinearSvmTrainer<InputType> base_type;

	LinearCSvmTrainer(double C, bool offset, bool unconstrained = false) 
	: AbstractLinearSvmTrainer<InputType>(C, offset, unconstrained){}

	/// \brief From INameable: return the class name.
	std::string name() const
	{ return "LinearCSvmTrainer"; }
	
	/// \brief sets the type of the multi-class svm used
	void setMcSvmType(McSvm type){
		m_McSvmType = type;
	}

	void train(LinearClassifier<InputType>& model, LabeledData<InputType, unsigned int> const& dataset)
	{
		std::size_t classes = numberOfClasses(dataset);
		if(classes == 2){
			trainBinary(model,dataset);
			return;
		}
		switch (m_McSvmType){
			case McSvm::WW:
				trainMc<QpMcLinearWW<InputType> >(model,dataset,classes);
			break;
			case McSvm::CS:
				trainMc<QpMcLinearCS<InputType> >(model,dataset,classes);
			break;
			case McSvm::LLW:
				trainMc<QpMcLinearLLW<InputType> >(model,dataset,classes);
			break;
			case McSvm::ATM:
				trainMc<QpMcLinearATM<InputType> >(model,dataset,classes);
			break;
			case McSvm::ATS:
				trainMc<QpMcLinearATS<InputType> >(model,dataset,classes);
			break;
			case McSvm::ADM:
				trainMc<QpMcLinearADM<InputType> >(model,dataset,classes);
			break;
			case McSvm::MMR:
				trainMc<QpMcLinearMMR<InputType> >(model,dataset,classes);
			break;
			case McSvm::ReinforcedSvm:
				trainMc<QpMcLinearReinforced<InputType> >(model,dataset,classes);
			break;
			case McSvm::OVA://OVA is a special case and implemented here
				trainOVA(model,dataset,classes);
			break;
		}
	}
private:
	McSvm m_McSvmType;

	void trainBinary(LinearClassifier<InputType>& model, LabeledData<InputType, unsigned int> const& dataset)
	{
		std::size_t dim = inputDimension(dataset);
		QpBoxLinear<InputType> solver(dataset, dim);
		solver.solve(
				base_type::C(),
				0.0,
				QpConfig::stoppingCondition(),
				&QpConfig::solutionProperties(),
				QpConfig::verbosity() > 0);
		
		if(!this->trainOffset()){
			RealMatrix w(1, dim, 0.0);
			row(w,0) = solver.solutionWeightVector();
			model.decisionFunction().setStructure(w);
			return;
		}
		
		double offset = 0;
		double stepSize = 0.1;
		double grad = solver.offsetGradient();
		while(stepSize > 0.1*QpConfig::stoppingCondition().minAccuracy){
			offset+= (grad < 0? -stepSize:stepSize);
			solver.setOffset(offset);
			solver.solve(
				base_type::C(),
				0.0,
				QpConfig::stoppingCondition(),
				&QpConfig::solutionProperties(),
				QpConfig::verbosity() > 0);
			double newGrad = solver.offsetGradient();
			if(newGrad == 0)
				break;
			if(newGrad*grad < 0)
				stepSize *= 0.5;
			else
				stepSize *= 1.6;
			grad = newGrad;
		}
		
		RealMatrix w(1, dim, 0.0);
		noalias(row(w,0)) = solver.solutionWeightVector();
		model.decisionFunction().setStructure(w,RealVector(1,offset));
		
	}
	template<class Solver>
	void trainMc(LinearClassifier<InputType>& model, LabeledData<InputType, unsigned int> const& dataset, std::size_t classes){
		std::size_t dim = inputDimension(dataset);

		Solver solver(dataset, dim, classes);
		RealMatrix w = solver.solve(random::globalRng, this->C(), this->stoppingCondition(), &this->solutionProperties(), this->verbosity() > 0);
		model.decisionFunction().setStructure(w);
	}
	
	void trainOVA(LinearClassifier<InputType>& model, const LabeledData<InputType, unsigned int>& dataset, std::size_t classes)
	{
		base_type::m_solutionproperties.type = QpNone;
		base_type::m_solutionproperties.accuracy = 0.0;
		base_type::m_solutionproperties.iterations = 0;
		base_type::m_solutionproperties.value = 0.0;
		base_type::m_solutionproperties.seconds = 0.0;

		std::size_t dim = inputDimension(dataset);
		RealMatrix w(classes, dim);
		for (unsigned int c=0; c<classes; c++)
		{
			LabeledData<InputType, unsigned int> bindata = oneVersusRestProblem(dataset, c);
			QpBoxLinear<InputType> solver(bindata, dim);
			QpSolutionProperties prop;
			solver.solve(this->C(), 0.0, base_type::m_stoppingcondition, &prop, base_type::m_verbosity > 0);
			noalias(row(w, c)) = solver.solutionWeightVector();
			base_type::m_solutionproperties.iterations += prop.iterations;
			base_type::m_solutionproperties.seconds += prop.seconds;
			base_type::m_solutionproperties.accuracy = std::max(base_type::solutionProperties().accuracy, prop.accuracy);
		}
		model.decisionFunction().setStructure(w);
	}
};


template <class InputType, class CacheType = float>
class SquaredHingeCSvmTrainer : public AbstractSvmTrainer<InputType, unsigned int>
{
public:
	typedef CacheType QpFloatType;

	typedef RegularizedKernelMatrix<InputType, QpFloatType> KernelMatrixType;
	typedef CachedMatrix< KernelMatrixType > CachedMatrixType;
	typedef PrecomputedMatrix< KernelMatrixType > PrecomputedMatrixType;

	typedef AbstractModel<InputType, RealVector> ModelType;
	typedef AbstractKernelFunction<InputType> KernelType;
	typedef AbstractSvmTrainer<InputType, unsigned int> base_type;

	//! Constructor
	//! \param  kernel         kernel function to use for training and prediction
	//! \param  C              regularization parameter - always the 'true' value of C, even when unconstrained is set
	//! \param  unconstrained  when a C-value is given via setParameter, should it be piped through the exp-function before using it in the solver??
	SquaredHingeCSvmTrainer(KernelType* kernel, double C, bool unconstrained = false)
	: base_type(kernel, C, unconstrained)
	{ }
	
	//! Constructor
	//! \param  kernel         kernel function to use for training and prediction
	//! \param  negativeC   regularization parameter of the negative class (label 0)
	//! \param  positiveC    regularization parameter of the positive class (label 1)
	//! \param  unconstrained  when a C-value is given via setParameter, should it be piped through the exp-function before using it in the solver?
	SquaredHingeCSvmTrainer(KernelType* kernel, double negativeC, double positiveC, bool unconstrained = false)
	: base_type(kernel,negativeC, positiveC, unconstrained)
	{ }

	/// \brief From INameable: return the class name.
	std::string name() const
	{ return "SquaredHingeCSvmTrainer"; }

	/// \brief Train the C-SVM.
	void train(KernelClassifier<InputType>& svm, LabeledData<InputType, unsigned int> const& dataset)
	{		
		svm.decisionFunction().setStructure(base_type::m_kernel,dataset.inputs(),this->m_trainOffset);
		
		RealVector diagonalModifier(dataset.numberOfElements(),0.5/base_type::m_regularizers(0));
		if(base_type::m_regularizers.size() != 1){
			for(std::size_t i = 0; i != diagonalModifier.size();++i){
				diagonalModifier(i) = 0.5/base_type::m_regularizers(dataset.element(i).label);
			}
		}
		
		KernelMatrixType km(*base_type::m_kernel, dataset.inputs(),diagonalModifier);
		if (QpConfig::precomputeKernel())
		{
			PrecomputedMatrixType matrix(&km);
			optimize(svm.decisionFunction(),matrix,diagonalModifier,dataset);
		}
		else
		{
			CachedMatrixType matrix(&km);
			optimize(svm.decisionFunction(),matrix,diagonalModifier,dataset);
		}
		base_type::m_accessCount = km.getAccessCount();
		if (base_type::sparsify()) svm.decisionFunction().sparsify();

	}

private:
	
	template<class Matrix>
	void optimize(KernelExpansion<InputType>& svm, Matrix& matrix,RealVector const& diagonalModifier, LabeledData<InputType, unsigned int> const& dataset){
		typedef CSVMProblem<Matrix> SVMProblemType;
		SVMProblemType svmProblem(matrix,dataset.labels(),1e100);
		if (this->m_trainOffset)
		{
			typedef SvmShrinkingProblem<SVMProblemType> ProblemType;
			ProblemType problem(svmProblem,base_type::m_shrinking);
			QpSolver< ProblemType > solver(problem);
			solver.solve(base_type::stoppingCondition(), &base_type::solutionProperties());
			column(svm.alpha(),0)= problem.getUnpermutedAlpha();
			//compute the bias
			double sum = 0.0;
			std::size_t freeVars = 0;
			for (std::size_t i=0; i < problem.dimensions(); i++)
			{
				if(problem.alpha(i) > problem.boxMin(i) && problem.alpha(i) < problem.boxMax(i)){
					sum += problem.gradient(i) - problem.alpha(i)*2*diagonalModifier(i);
					freeVars++;
				}
			}
			if (freeVars > 0)
				svm.offset(0) =  sum / freeVars;		//stabilized (averaged) exact value
			else
				svm.offset(0) = 0;
		}
		else
		{
			typedef BoxConstrainedShrinkingProblem<SVMProblemType> ProblemType;
			ProblemType problem(svmProblem,base_type::m_shrinking);
			QpSolver< ProblemType > solver(problem);
			solver.solve(base_type::stoppingCondition(), &base_type::solutionProperties());
			column(svm.alpha(),0) = problem.getUnpermutedAlpha();
			
		}
	}
};


template <class InputType>
class SquaredHingeLinearCSvmTrainer : public AbstractLinearSvmTrainer<InputType>
{
public:
	typedef AbstractLinearSvmTrainer<InputType> base_type;

	SquaredHingeLinearCSvmTrainer(double C, bool unconstrained = false) 
	: AbstractLinearSvmTrainer<InputType>(C, false, unconstrained){}

	/// \brief From INameable: return the class name.
	std::string name() const
	{ return "SquaredHingeLinearCSvmTrainer"; }

	void train(LinearClassifier<InputType>& model, LabeledData<InputType, unsigned int> const& dataset)
	{
		std::size_t dim = inputDimension(dataset);
		QpBoxLinear<InputType> solver(dataset, dim);
		RealMatrix w(1, dim, 0.0);
		solver.solve(
				1e100,
				1.0 / base_type::C(),
				QpConfig::stoppingCondition(),
				&QpConfig::solutionProperties(),
				QpConfig::verbosity() > 0);
		row(w,0) = solver.solutionWeightVector();
		model.decisionFunction().setStructure(w);
	}
};


}
#endif