adaboost.cpp

AdaBoost is an efficient tool in machine learning. It can combine a series of weak learners into a s

/****************************************************************************
 NJU Magic.  Copyright (c) 2007.  All Rights Reserved.            
  
--------------------------------------------------------------------
Permission to use, copy, or modify this software and its documentation
for educational and research purposes only and without fee is hereby
granted, provided that this copyright notice appear on all copies and
supporting documentation.  For any other uses of this software, in
original or modified form, including but not limited to distribution
in whole or in part, specific prior permission must be obtained from
NJU Magic and the authors.  These programs shall not be used, rewritten,
or adapted as the basis of a commercial software or hardware product
without first obtaining appropriate licenses from NJU Magic. NJU Magic 
makes no representations about the suitability of this software for any
purpose.  It is provided "as is" without express or implied warranty.

---------------------------------------------------------------------

 File: AdaBoost.cpp                                                                         
 Authors: Yao Wei  
 Date Created    : 2007-8-11
                                          
****************************************************************************/

#include "AdaBoost.h"

//Construction
Boosting::Boosting():leaner(0),	max_iter(0),n_input(0),
	n_samples(0), submat(0)
{
	
}

// Destruction
Boosting::~Boosting()
{
	Free();
}

void Boosting::Preprocess(const Mat &train, const Vec &responses)
{
	if (train.rows() != responses.length())
	{
		mexPrintf("ERROR: # train sample != # lable\n");
		exit(-1);
	}
	
	n_samples = train.rows();
	n_input = train.cols();
	
	GetClassInfo(responses);
	UnrollData(train,responses);
}

void Boosting::GetClassInfo(const Vec &responses)
{
	int i = 0, j = 0, k = 0;
	double a[100];				//store class label
	int b[100] = {0};			//each label's count
	double temp;
	int count = 0;				//#class
	
	while(i < responses.length())
	{
		temp = responses[i];	
		k = 0;
		while(k < j && temp != a[k])
		{
			k++;
		}
		if(k == j)				//a new class comes
		{
			count ++;
			b[j] ++;			//correspond count+1
			a[j] = temp;
			j ++;
		}
		else					//the class has been existed
		{
			b[k] ++;			//correspond count+1
		}
		i ++;
	}

	//get # class
	info.count = count;
	//get count of each class label
	info.eachcount = new int[count];
	for(i = 0; i < count; i ++)
		info.eachcount[i] = b[i];
	//get concrete class label
	info.eachlabel = new double[count];
	for(i = 0; i < count; i ++)
		info.eachlabel[i] = a[i];

	/*
	mexPrintf("Training data consists of %d classes.\n", count);
	mexPrintf("They are ");
	for(i = 0; i< count; i++)
		mexPrintf("%g(#%d) ", a[i], b[i]);
	mexPrintf("repectively.\n");
	*/
}

void Boosting::UnrollData(const Mat &train,const Vec &responses)
{
	int i, j;
	int count = info.count;
	
	//1 vs 1, get sub-matrix(s) according to each class
	//that is to say, each sample in a submat shares the same label
	//if # class if count, we will get count submat(s) finally
	submat = new Mat[count];
	for(i=0; i<count; i++)
		submat[i].Set(info.eachcount[i],n_input);
	
	for(i=0; i<n_samples; i++)
	{
		for(j=0; j<count; j++)
		{
			if(responses[i] == info.eachlabel[j])
			{
				submat[j].AddRow(train[i]);
				break;
			}
		}
	}
}

void Boosting::DoTrain()
{
	int count = info.count;			//#	class
	leaner = new WeakLearner[count*(count-1)/2];
	int k=0;
	//1 vs 1, do count*(count-1)/2 steps of binary training
	/*
	mexPrintf("\n------Training Begin------");
	mexPrintf("\nThe whole training process contains %d steps of binary training, "
		"each training iterarion is %d.", 
		count*(count-1)/2, max_iter);
	*/

	for(int i=0; i<count; i++)
	{
		for(int j=i+1; j<count; j++)
		{
			//mexPrintf("\nstep %3d: %g vs %g\n", k+1, 
			//	info.eachlabel[i], info.eachlabel[j]);
			leaner[k] = TrainOneVsOne(submat[i], submat[j],
				info.eachlabel[i], info.eachlabel[j]);
			k++;
		}
		submat[i].Free();
	}
	//mexPrintf("\n------Training OK------\n\n");
	delete []submat; submat = 0;
	
}

WeakLearner Boosting::
TrainOneVsOne(const Mat&mat1, const Mat&mat2,double label1,double label2)
{
	int i, t;
	int n = mat1.rows() + mat2.rows();
	int D = mat1.cols();
	double* _weights = new double[n];
	// Initialize weights
	for(i = 0; i< n; i++)
		_weights[i] = 1.0 / n;
	int* hClassification = new int[n];

	WeakLearner wlearner;
	wlearner.Init(max_iter, label1, label2);
	DecisionStump stump;
	
	// Perform the learning
	for(t = 0; t < max_iter; t ++)
	{
		//if(t % 10 == 0)	mexPrintf("*...\r");
		
		// Create the weak learner and train it
		stump.RoundTrain(mat1, mat2, _weights);
		
		// Compute the classifications and training error
		double epsilon = 0.0;
		for(i = 0; i < mat1.rows(); i ++)
		{
			hClassification[i] = stump.Classify(mat1[i], D);
			epsilon += (hClassification[i] == +1) ? 0: _weights[i];
		}
		for(i = 0; i< mat2.rows(); i++)
		{
			hClassification[i+mat1.rows()] = stump.Classify(mat2[i], D);
			epsilon += 
				(hClassification[i+mat1.rows()] == -1) ? 0: _weights[i+mat1.rows()];
		}

		// Check stopping condition
		//if(epsilon >= 0.5)	break;
		
		// Calculate alpha
		double alpha  = 0.5 *log((1 - epsilon) / epsilon);

		// Update the weights
		double weightsSum = 0.0;
		for (i = 0; i < n; i++)
		{
			_weights[i] *= 
				exp(-alpha * (i < mat1.rows() ? +1 : -1) * hClassification[i]);
			weightsSum += _weights[i];
		}

		// Normalize
		for (i = 0; i < n; i++)
			_weights[i] /= weightsSum;

		// Store the weak learner and alpha value
		wlearner.AddStump(stump, alpha);
	}
	delete []_weights;
	delete []hClassification;

	return wlearner;
}

double Boosting::Vote(double *feature,int n_input)
{
	if(this->n_input != n_input)
	{
		mexPrintf("The feature of test file doesn't match to the model file!\n");
		exit(-1);
	}	
	
	int count = info.count;
	// dicision[] which store the count of possible 
	// class according to the weaklearner  
	int *dicision = new int[count];
	memset(dicision, 0, count*sizeof(int));
	int i,j;
	double tempresult;
	// Count the predict result of each weaklearner
	for(i=0; i<count*(count-1)/2; i++)
	{	
		tempresult = leaner[i].Classify(feature, n_input);
		for(j=0; j<count; j++)
		{
			if(info.eachlabel[j] == tempresult)
			{
				dicision[j]++;
				break;
			}
		}
	}

	// majority vote
	double maxcount = dicision[0];
	double label = info.eachlabel[0];
	for(i=1; i<count; i++)
	{
		if(maxcount < dicision[i])
		{
			label = info.eachlabel[i];
			maxcount = dicision[i];
		}
	}
	
	return label;			
}



void Boosting::Free()
{
	for(int i = 0; i < info.count*(info.count-1)/2; i++)
		leaner[i].Free();
	delete []leaner; leaner = 0;
	delete []info.eachcount; info.eachcount = 0;
	delete []info.eachlabel; info.eachlabel = 0;
}

/************************************************************************/
/* For Train                                                            */
/************************************************************************/
// Input Arguments
#define TRAINDATA prhs[0]
#define TRAINLABEL prhs[1]
#define TRAINSTEP prhs[2]

// Out Arguments
#define TRAINMODEL plhs[0]


void Boosting::mexTrain(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
	if(nrhs <2 || nrhs >3)
	{
		mexPrintf("Usage: model = boosttrain(train_data, train_label, iteration);\n\n");
		return;
	}
	int i, j;
	/* get the size of train_data */
	int m, n;
	const int *m_n = mxGetDimensions(TRAINDATA);
	m = m_n[0];
	n = m_n[1];
	
	/* convert Matlab data & label to C++ */
	Mat train_data(m, n);
	double * data = (double*)mxGetData(TRAINDATA);
	for(i = 0; i < m; i++)
		for(j = 0; j < n ; j++)
			train_data[i][j] = data[j*m + i];
	
	Vec train_label(m);
	double * label = (double*)mxGetData(TRAINLABEL);
	for(i = 0; i < m ; i++)
		train_label[i] = label[i];

	if(nrhs == 2)
		max_iter = 100;
	else
	{
		double *iter = (double*)mxGetData(TRAINSTEP);
		max_iter = (int)iter[0];
	}
	
	/* Run Boost Train */
	Preprocess(train_data, train_label);
	train_data.Free();
	train_label.Free();
	DoTrain();

	/* convert C++ BoostModel to Matlab */
	const char *field_names[]={
		"n_size",
		"n_step",
		"nr_class",
		"totalWLs",
		"Labels",
		"WLs"
	};
	int num_of_fields = 6;

	mxArray **rhs;
	double *ptr;
	rhs = (mxArray**)mxMalloc(sizeof(mxArray*)*num_of_fields);
	int out_id = 0;
	
	//n_size
	rhs[out_id] = mxCreateDoubleMatrix(2,1,mxREAL);
	ptr = mxGetPr(rhs[out_id]);
	ptr[0] = m;
	ptr[1] = n;
	out_id++;

	//n_step
	rhs[out_id] = mxCreateDoubleMatrix(1,1,mxREAL);
	ptr = mxGetPr(rhs[out_id]);
	ptr[0] = max_iter;
	out_id ++;
	
	//nr_class
	rhs[out_id] = mxCreateDoubleMatrix(1,1,mxREAL);
	ptr = mxGetPr(rhs[out_id]);
	ptr[0] = info.count;
	out_id++;

	//totalWLs
	rhs[out_id] = mxCreateDoubleMatrix(1,1,mxREAL);
	ptr = mxGetPr(rhs[out_id]);
	ptr[0] = info.count*(info.count-1)/2;
	out_id++;

	//Labels
	rhs[out_id] = mxCreateDoubleMatrix(info.count,1,mxREAL);
	ptr = mxGetPr(rhs[out_id]);
	for(i = 0; i < info.count; i++)
		ptr[i] = info.eachlabel[i];
	out_id++;

	//WLs
	rhs[out_id] = mxCreateDoubleMatrix(4*max_iter,info.count*(info.count-1)/2,
		mxREAL);
	ptr = mxGetPr(rhs[out_id]);
	for(i = 0; i < info.count*(info.count-1)/2; i++)
	{
		for(j = 0; j < max_iter; j++)
		{
			ptr[4*max_iter*i+4*j+0] = leaner[i].alpha[j];
			ptr[4*max_iter*i+4*j+1] = leaner[i].decisionstump[j]._d;
			ptr[4*max_iter*i+4*j+2] = leaner[i].decisionstump[j]._sign;
			ptr[4*max_iter*i+4*j+3] = leaner[i].decisionstump[j]._threshold;
		}
		
	}

	/* Create a struct matrix contains NUM_OF_RETURN_FIELD fields */
	mxArray *return_model = mxCreateStructMatrix(1, 1, num_of_fields, field_names);

	/* Fill struct matrix with input arguments */
	for(i = 0; i < num_of_fields; i++)
		mxSetField(return_model,0,field_names[i],mxDuplicateArray(rhs[i]));
	
	/* return */
	TRAINMODEL = return_model;
	
	mxFree(rhs);
}

/************************************************************************/
/* For Predict                                                          */
/************************************************************************/
// Input Arguments
#define MODEL prhs[0]
#define PREDICTDATA prhs[1]
#define TESTLABEL prhs[2]

// Out Arguments
#define PREDICTLABEL plhs[0]
#define ACCURACY plhs[1]

void Boosting::mexPredict(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
	if(nrhs < 2 || nrhs >3)
	{
		mexPrintf("Usage: predict_label = boostpredict(model, predict_data);\n");
		mexPrintf("   or  [predict_label, accuracy] = boostpredict(model, test_data, test_label);\n\n");
		return;
	}

	int i, j, k, l;
	
	/* get the size of predict_data */
	int m, n;
	const int *m_n = mxGetDimensions(PREDICTDATA);
	m = m_n[0];
	n = m_n[1];

	/* convert Matlab-Model to C++ */
	int id = 0;
	double *ptr;
	mxArray **rhs;
	int num_of_fields;
	num_of_fields = mxGetNumberOfFields(MODEL);
	rhs = (mxArray **) mxMalloc(sizeof(mxArray *)*num_of_fields);
	for(i=0;i<num_of_fields;i++)
		rhs[i] = mxGetFieldByNumber(MODEL, 0, i);

	//n_size
	ptr = mxGetPr(rhs[id]);
	n_samples = (int)ptr[0];
	n_input  = (int)ptr[1];
	id++;

	//n_step
	ptr = mxGetPr(rhs[id]);
	max_iter = (int)ptr[0];
	id++;

	//nr_class
	ptr = mxGetPr(rhs[id]);
	info.count = (int)ptr[0];
	id++;

	//totalWLs
	ptr = mxGetPr(rhs[id]);
	int n_WLs = (int)ptr[0];
	id++;

	//Labels
	ptr = mxGetPr(rhs[id]);
	info.eachlabel = new double[info.count];
	for(i = 0; i < info.count; i++)
		info.eachlabel[i] = ptr[i];
	id++;
	
	//WLs
	ptr = mxGetPr(rhs[id]);
	leaner = new WeakLearner[n_WLs];
	for(i = 0; i < n_WLs; i++)
	{
		leaner[i].alpha = new double[max_iter];
		memset(leaner[i].alpha, 0, sizeof(double)*max_iter);
		leaner[i].decisionstump = new DecisionStump[max_iter];
	}
	i = 0;
	for(k = 0; k < info.count; k++)
	{
		for(l = k+1; l < info.count; l++)
		{
			for(j = 0; j < max_iter; j++)
			{
				leaner[i].alpha[j] = ptr[4*max_iter*i+4*j+0];
				leaner[i].decisionstump[j]._d = (int)ptr[4*max_iter*i+4*j+1];
				leaner[i].decisionstump[j]._sign = (int)ptr[4*max_iter*i+4*j+2];
				leaner[i].decisionstump[j]._threshold = ptr[4*max_iter*i+4*j+3];
				leaner[i].maxnum = max_iter;
				leaner[i]._iterator = max_iter;
				leaner[i].positive_label = info.eachlabel[k];
				leaner[i].negative_label = info.eachlabel[l];
			}

			i++;
		}
	}	
		
	/* Allocate memory for output result */
	PREDICTLABEL = mxCreateDoubleMatrix(m, 1, mxREAL);
	ptr = mxGetPr(PREDICTLABEL);
	
	/* one sample features used to predict */
	double * feature = new double[n];

	/* predict each sample */
	double * data = (double*)mxGetData(PREDICTDATA);
	for(i = 0; i < m; i++)
	{
		for(j = 0; j < n ; j++)
			feature[j] = data[j*m + i];
		//predict its label
		ptr[i] = Vote(feature, n);
	}
	/* Computer the classify accuracy*/
	if(nrhs == 3)
	{
		int sum = 0;
		ptr = (double*)mxGetData(TESTLABEL);
		data = (double*)mxGetData(PREDICTLABEL);
		for(i = 0; i < m; i ++)
		{
			if(ptr[i] == data[i])
				sum++;
		}
		ACCURACY = mxCreateDoubleMatrix(1, 1, mxREAL);
		ptr = mxGetPr(ACCURACY);
		ptr[0] = (double)sum / m;
	}

	delete []feature;
	mxFree(rhs);

}