added documentation, restart git usage, start implementing kernels

18 files changed, 1878 insertions, 423 deletions

diff --git a/src/crossval.c b/src/crossval.c
index 9a3c1cc..10e3051 100644
--- a/src/crossval.c
+++ b/src/crossval.c
@@ -1,7 +1,40 @@
+/**
+ * @file crossval.c
+ * @author Gertjan van den Burg 
+ * @date January 7, 2014
+ * @brief Functions for cross validation
+ *
+ * @details
+ * This file contains functions for performing cross validation. The funtion
+ * msvmmaj_make_cv_split() creates a cross validation vector for non-stratified
+ * cross validation. The function msvmmaj_get_tt_split() creates a train and 
+ * test dataset from a given dataset and a pre-determined CV partition vector.
+ * See individual function documentation for details.
+ *
+ */
+
 #include "crossval.h"
-#include "matrix.h"
-#include "MSVMMaj.h"
+#include "msvmmaj.h"
+#include "msvmmaj_matrix.h"
 
+/**
+ * @brief Create a cross validation split vector
+ *
+ * @details
+ * A pre-allocated vector of length N is created which can be used to define
+ * cross validation splits. The folds are contain between 
+ * @f$ \lfloor N / folds \rfloor @f$ and @f$ \lceil N / folds \rceil @f$ 
+ * instances. An instance is mapped to a partition randomly until all folds 
+ * contain @f$ N \% folds @f$ instances. The zero fold then contains 
+ * @f$ N / folds + N \% folds @f$ instances. These remaining @f$ N \% folds @f$
+ * instances are then distributed over the first @f$ N \% folds @f$ folds. 
+ *
+ * @param[in] 		N 	number of instances
+ * @param[in] 		folds 	number of folds
+ * @param[in,out] 	cv_idx 	array of size N which contains the fold index
+ * 				for each observation on exit	
+ *
+ */
 void msvmmaj_make_cv_split(long N, long folds, long *cv_idx)
 {
 	long i, j, idx;
@@ -30,6 +63,26 @@ void msvmmaj_make_cv_split(long N, long folds, long *cv_idx)
 	}
 }
 
+
+/**
+ * @brief Create train and test datasets for a CV split
+ *
+ * @details
+ * Given a MajData structure for the full dataset, a previously created
+ * cross validation split vector and a fold index, a training and test dataset
+ * are created. 
+ *
+ * @param[in] 		full_data 	a MajData structure for the entire 
+ * 					dataset
+ * @param[in,out] 	train_data 	an initialized MajData structure which
+ * 					on exit contains the training dataset
+ * @param[in,out] 	test_data 	an initialized MajData structure which
+ * 					on exit contains the test dataset
+ * @param[in] 		cv_idx 		a vector of cv partitions created by
+ * 					msvmmaj_make_cv_split()
+ * @param[in] 		fold_idx 	index of the fold which becomes the 
+ * 					test dataset
+ */
 void msvmmaj_get_tt_split(struct MajData *full_data, struct MajData *train_data,
 		struct MajData *test_data, long *cv_idx, long fold_idx)
 {
@@ -67,13 +120,15 @@ void msvmmaj_get_tt_split(struct MajData *full_data, struct MajData *train_data,
 			test_data->y[k] = full_data->y[i];
 			for (j=0; j<m+1; j++)
 				matrix_set(test_data->Z, m+1, k, j, 
-						matrix_get(full_data->Z, m+1, i, j));
+						matrix_get(full_data->Z, m+1, 
+							i, j));
 			k++;
 		} else {
 			train_data->y[l] = full_data->y[i];
 			for (j=0; j<m+1; j++)
 				matrix_set(train_data->Z, m+1, l, j,
-						matrix_get(full_data->Z, m+1, i, j));
+						matrix_get(full_data->Z, m+1, 
+							i, j));
 			l++;
 		}
 	}
diff --git a/src/kernel.c b/src/kernel.c
deleted file mode 100644
index ee64871..0000000
--- a/src/kernel.c
+++ /dev/null
@@ -1,85 +0,0 @@
-/**
- * @file kernel.c
- * @author Gertjan van den Burg (burg@ese.eur.nl)
- * @date October 18, 2013
- * @brief Defines main functions for use of kernels in MSVMMaj.
- *
- * @details
- * Functions for constructing different kernels using user-supplied 
- * parameters. Also contains the functions for decomposing the
- * kernel matrix using several decomposition methods.
- *
- */
-#include <math.h>
-
-#include "kernel.h"
-
-void msvmmaj_make_kernel(struct MajModel *model, struct MajData *data)
-{
-	switch (model->kerneltype) {
-		case K_LINEAR:
-			break;
-		case K_POLY:
-			msvmmaj_make_kernel_poly(model, data);
-			break;
-		case K_RBF:
-			msvmmaj_make_kernel_rbf(model, data);
-			break;
-		case K_SIGMOID:
-			msvmmaj_make_kernel_sigmoid(model, data);
-			break;
-	}
-}
-
-void msvmmaj_make_kernel_rbf(struct MajModel *model, struct MajData *data)
-{
-	long i, j;
-	long n = model->n;
-	double value;
-	double *x1, *x2;
-	double *K = Calloc(double, n*(n+1));
-
-	for (i=0; i<n; i++) {
-		for (j=0; j<n; j++) {
-			x1 = &data->Z[i*(data->m+1)+1];
-			x2 = &data->Z[j*(data->m+1)+1];
-			value = msvmmaj_compute_rbf(x1, x2, model->kernelparam, n);
-			matrix_set(K, n+1, i, j+1, value);
-		}
-		matrix_set(K, n+1, i, 0, 1.0);
-	}
-
-	free(data->Z);
-	data->Z = K;
-	data->m = n;
-	model->m = n;
-}
-
-/**
- * Implements k(x, z) = exp( -gamma * || x - z ||^2)
- */
-double msvmmaj_compute_rbf(double *x1, double *x2, double *kernelparam, long n)
-{
-	long i;
-	double value = 0.0;
-
-	for (i=0; i<n; i++) 
-		value += (x1[i] - x2[i]) * (x1[i] - x2[i]);
-	value *= -kernelparam[0];
-	return exp(value);
-}
-
-/**
- * Implements k(x, z) = (gamma * <x, z> + c)^degree
- */
-double msvmmaj_compute_poly(double *x1, double *x2, double *kernelparam, long n)
-{
-	long i;
-	double value = 0.0;
-	for (i=0; i<n; i++)
-		value += x1[i]*x2[i];
-	value *= kernelparam[0];
-	value += kernelparam[1];
-	for (i=1; i<(int kernelparam[2]); i++)
-		value *= value;
-	:w
diff --git a/src/libMSVMMaj.c b/src/libMSVMMaj.c
index 9544830..a0bef97 100644
--- a/src/libMSVMMaj.c
+++ b/src/libMSVMMaj.c
@@ -1,6 +1,6 @@
 /**
  * @file libMSVMMaj.c
- * @author Gertjan van den Burg (burg@ese.eur.nl)
+ * @author Gertjan van den Burg
  * @date August 8, 2013
  * @brief Main functions for the MSVMMaj algorithm
  *
@@ -16,24 +16,23 @@
 #include <math.h>
 
 #include "libMSVMMaj.h"
-#include "MSVMMaj.h"
-#include "matrix.h"
+#include "msvmmaj.h"
+#include "msvmmaj_matrix.h"
 
 inline double rnd() { return (double) rand()/0x7FFFFFFF; }
 
 /**
- * @name msvmmaj_simplex_gen
  * @brief Generate matrix of simplex vertex coordinates
- * @ingroup libMSVMMaj
  * 
+ * @details
  * Generate the simplex matrix. Each row of the created 
  * matrix contains the coordinate vector of a single 
  * vertex of the K-simplex in K-1 dimensions. The simplex
  * generated is a special simplex with edges of length 1.
  * The simplex matrix U must already have been allocated.
  *
- * @param [in] 		K 	number of classes
- * @param [in,out] 	U 	simplex matrix of size K * (K-1)
+ * @param[in] 		K 	number of classes
+ * @param[in,out] 	U 	simplex matrix of size K * (K-1)
  */
 void msvmmaj_simplex_gen(long K, double *U)
 {
@@ -51,10 +50,18 @@ void msvmmaj_simplex_gen(long K, double *U)
 	}
 }
 
-/*!
-	Generate the category matrix R. The category matrix has 1's everywhere
-	except at the column corresponding to the label of instance i.
-*/
+/**
+ * @brief Generate the category matrix
+ *
+ * @details
+ * Generate the category matrix R. The category matrix has 1's everywhere
+ * except at the column corresponding to the label of instance i, there the
+ * element is 0.
+ *
+ * @param[in,out] 	model 		corresponding MajModel
+ * @param[in] 		dataset 	corresponding MajData
+ *
+ */
 void msvmmaj_category_matrix(struct MajModel *model, struct MajData *dataset)
 {
 	long i, j;
@@ -70,8 +77,19 @@ void msvmmaj_category_matrix(struct MajModel *model, struct MajData *dataset)
 	}
 }
 
-/*!
- * Simplex diff
+/**
+ * @brief Generate the simplex difference matrix
+ *
+ * @details
+ * The simplex difference matrix is a 3D matrix which is constructed
+ * as follows. For each instance i, the difference vectors between the row of 
+ * the simplex matrix corresponding to the class label of instance i and the
+ * other rows of the simplex matrix are calculated. These difference vectors 
+ * are stored in a matrix, which is one horizontal slice of the 3D matrix.
+ *
+ * @param[in,out] 	model 	the corresponding MajModel
+ * @param[in] 		data 	the corresponding MajData
+ *
  */
 void msvmmaj_simplex_diff(struct MajModel *model, struct MajData *data)
 {
@@ -92,13 +110,22 @@ void msvmmaj_simplex_diff(struct MajModel *model, struct MajData *data)
 	}
 }
 
-/*!
-	Calculate the errors Q based on the current value of V.
-	It is assumed that the memory for Q has already been allocated. 
-	In addition, the matrix ZV is calculated here. It is assigned to a 
-	pre-allocated block of memory, since it would be inefficient to keep
-	reassigning this block at every iteration.
-*/
+/**
+ * @brief Calculate the scalar errors
+ * 
+ * @details
+ * Calculate the scalar errors q based on the current estimate of V, and
+ * store these in Q. It is assumed that the memory for Q has already been 
+ * allocated. In addition, the matrix ZV is calculated here. It is assigned 
+ * to a pre-allocated block of memory, which is passed to this function.
+ *
+ * @param[in,out] 	model 	the corresponding MajModel
+ * @param[in] 		data 	the corresponding MajData
+ * @param[in,out] 	ZV 	a pointer to a memory block for ZV. On exit
+ * 				this block is updated with the new ZV matrix
+ * 				calculated with MajModel::V.
+ *
+ */
 void msvmmaj_calculate_errors(struct MajModel *model, struct MajData *data, double *ZV)
 {
 	long i, j, k;
@@ -136,9 +163,23 @@ void msvmmaj_calculate_errors(struct MajModel *model, struct MajData *data, doub
 	}
 }	
 
-/*!
-	Calculate the Huber hinge errors for each error in the matrix Q.
-*/
+/**
+ * @brief Calculate the Huber hinge errors
+ *
+ * @details
+ * For each of the scalar errors in Q the Huber hinge errors are 
+ * calculated. The Huber hinge is here defined as 
+ * @f[
+ * 	h(q) = 
+ * 		\begin{dcases}
+ * 			1 - q - \frac{\kappa + 1}{2} & \text{if } q \leq -\kappa \\
+ * 			\frac{1}{2(\kappa + 1)} ( 1 - q)^2 & \text{if } q \in (-\kappa, 1] \\
+ * 			0 & \text{if } q > 1
+ * 		\end{dcases}
+ * @f]
+ *
+ * @param[in,out] model 	the corresponding MajModel
+ */
 void msvmmaj_calculate_huber(struct MajModel *model) 
 {
 	long i, j;
@@ -159,10 +200,9 @@ void msvmmaj_calculate_huber(struct MajModel *model)
 }
 
 /**
- * @name msvmmaj_seed_model_V
  * @brief seed the matrix V from an existing model or using rand
- * @ingroup libMSVMMaj
  *
+ * @details
  * The matrix V must be seeded before the main_loop() can start. 
  * This can be done by either seeding it with random numbers or 
  * using the solution from a previous model on the same dataset
@@ -170,8 +210,8 @@ void msvmmaj_calculate_huber(struct MajModel *model)
  * significant improvement in the number of iterations necessary
  * because the seeded model V is closer to the optimal V.
  *
- * @param [in] 	from_model 	model from which to copy V
- * @param [in,out] to_model 	model to which V will be copied
+ * @param[in] 		from_model 	MajModel from which to copy V
+ * @param[in,out] 	to_model 	MajModel to which V will be copied
  */
 void msvmmaj_seed_model_V(struct MajModel *from_model, struct MajModel *to_model)
 {
@@ -193,10 +233,17 @@ void msvmmaj_seed_model_V(struct MajModel *from_model, struct MajModel *to_model
 	}
 }
 
-/*!
- * Step doubling
+/**
+ * @brief Use step doubling
+ *
+ * @details
+ * Step doubling can be used to speed up the Majorization algorithm. Instead 
+ * of using the value at the minimimum of the majorization function, the value
+ * ``opposite'' the majorization point is used. This can essentially cut the
+ * number of iterations necessary to reach the minimum in half.
+ *
+ * @param[in] 	model	MajModel containing the augmented parameters
  */
-
 void msvmmaj_step_doubling(struct MajModel *model)
 {
 	long i, j;
@@ -207,15 +254,33 @@ void msvmmaj_step_doubling(struct MajModel *model)
 	for (i=0; i<m+1; i++) {
 		for (j=0; j<K-1; j++) {
 			matrix_mul(model->V, K-1, i, j, 2.0);
-			matrix_add(model->V, K-1, i, j, -matrix_get(model->Vbar, K-1, i, j));
+			matrix_add(model->V, K-1, i, j, 
+					-matrix_get(model->Vbar, K-1, i, j));
 		}
 	}
 }
 
-/*!
- * initialize_weights
+/**
+ * @brief Initialize instance weights
+ *
+ * @details
+ * Instance weights can for instance be used to add additional weights to 
+ * instances of certain classes. Two default weight possibilities are
+ * implemented here. The first is unit weights, where each instance gets 
+ * weight 1. 
+ *
+ * The second are group size correction weights, which are calculated as
+ * @f[
+ * 	\rho_i = \frac{n}{Kn_k} ,
+ * @f]
+ * where @f$ n_k @f$ is the number of instances in group @f$ k @f$ and
+ * @f$ y_i = k @f$.
+ *
+ * @param[in] 		data 	MajData with the dataset
+ * @param[in,out] 	model 	MajModel with the weight specification. On 
+ * 				exit MajModel::rho contains the instance 
+ * 				weights.
  */
-
 void msvmmaj_initialize_weights(struct MajData *data, struct MajModel *model)
 {
 	long *groups;
diff --git a/src/matrix.c b/src/matrix.c
deleted file mode 100644
index 8803e8b..0000000
--- a/src/matrix.c
+++ /dev/null
@@ -1,77 +0,0 @@
-/**
- * @file matrix.c
- * @author Gertjan van den Burg (burg@ese.eur.nl)
- * @date August 8, 2013
- * @brief Functions facilitating matrix access
- *
- * @details
- * The functions contained in this file are used when
- * accessing or writing to matrices. Seperate functions
- * exist of adding and multiplying existing matrix
- * elements, to ensure this is done in place.
- *
- */
-
-#include "matrix.h"
-#include "util.h"
-
-/**
- * @name matrix_set
- * @brief Set element of matrix 
- * @ingroup matrix
- *
- * Row-Major order is used to set a matrix element. Since matrices
- * of type double are most common in MSVMMaj, this function only
- * deals with that type.
- *
- * @param [in] 	M 	matrix to set element of
- * @param [in]  cols 	number of columns of M
- * @param [in] 	i 	row index of element to write to
- * @param [in] 	j 	column index of element to write to
- * @param [out] val 	value to write to specified element of M
- */
-void matrix_set(double *M, long cols, long i, long j, double val)
-{
-	M[i*cols+j] = val;
-}
-
-double matrix_get(double *M, long cols, long i, long j)
-{
-	return M[i*cols+j];
-}
-
-void matrix_add(double *M, long cols, long i, long j, double val)
-{
-	M[i*cols+j] += val;
-}
-
-void matrix_mul(double *M, long cols, long i, long j, double val)
-{
-	M[i*cols+j] *= val;
-}
-
-void matrix3_set(double *M, long N2, long N3, long i, long j, 
-		long k, double val)
-{
-	M[k+N3*(j+N2*i)] = val;
-}
-
-double matrix3_get(double *M, long N2, long N3, long i, long j,
-		long k)
-{
-	return M[k+N3*(j+N2*i)];
-}
-
-
-void print_matrix(double *M, long rows, long cols)
-{
-	long i, j;
-
-	for (i=0; i<rows; i++) {
-		for (j=0; j<cols; j++)
-			note("%8.8f ", matrix_get(M, cols, i, j));
-		note("\n");
-	}
-	note("\n");
-}
-
diff --git a/src/msvmmaj_init.c b/src/msvmmaj_init.c
new file mode 100644
index 0000000..14278f9
--- /dev/null
+++ b/src/msvmmaj_init.c
@@ -0,0 +1,64 @@
+/**
+ * @file msvmmaj_init.c
+ * @author Gertjan van den Burg
+ * @date January 7, 2014
+ * @brief Functions for initializing model and data structures
+ *
+ * @details 
+ * This file contains functions for initializing a MajModel instance
+ * and a MajData instance. In addition, default values for these 
+ * structures are defined here (and only here).
+ *
+ */
+
+#include <math.h>
+
+#include "msvmmaj.h"
+#include "msvmmaj_init.h"
+
+/**
+ * @brief Initialize a MajModel structure
+ *
+ * @details
+ * A MajModel structure is initialized and the default value for the
+ * parameters are set. A pointer to the initialized model is returned.
+ *
+ * @returns 	initialized MajModel
+ */
+struct MajModel *msvmmaj_init_model()
+{
+	struct MajModel *model = Malloc(struct MajModel, 1);
+	
+	// set default values
+	model->p = 1.0;
+	model->lambda = pow(2, -8.0);
+	model->epsilon = 1e-6;
+	model->kappa = 0.0;
+	model->weight_idx = 1;
+	model->kerneltype = K_LINEAR;
+	model->use_cholesky = false;
+
+	return model;
+}
+
+/**
+ * @brief Initialize a MajData structure
+ * 
+ * @details
+ * A MajData structure is initialized and default values are set. 
+ * A pointer to the initialized data is returned.
+ *
+ * @returns 	initialized MajData
+ *
+ */
+struct MajData *msvmmaj_init_data()
+{
+	struct MajData *data = Malloc(struct MajData, 1);
+
+	// set default values
+	data->kerneltype = K_LINEAR;
+	data->use_cholesky = false;
+
+	return data;
+}	
+
diff --git a/src/msvmmaj_kernel.c b/src/msvmmaj_kernel.c
new file mode 100644
index 0000000..6238fc1
--- /dev/null
+++ b/src/msvmmaj_kernel.c
@@ -0,0 +1,195 @@
+/**
+ * @file msvmmaj_kernel.c
+ * @author Gertjan van den Burg
+ * @date October 18, 2013
+ * @brief Defines main functions for use of kernels in MSVMMaj.
+ *
+ * @details
+ * Functions for constructing different kernels using user-supplied 
+ * parameters. Also contains the functions for decomposing the
+ * kernel matrix using several decomposition methods.
+ *
+ */
+#include <math.h>
+
+#include "msvmmaj.h"
+#include "msvmmaj_kernel.h"
+#include "msvmmaj_lapack.h"
+#include "msvmmaj_matrix.h"
+#include "util.h"
+
+/**
+ * @brief Create the kernel matrix
+ *
+ * Create a kernel matrix based on the specified kerneltype. Kernel parameters 
+ * are assumed to be specified in the model.
+ *
+ * @param[in] 	model 	MajModel specifying the parameters
+ * @param[in] 	data 	MajData specifying the data.
+ *
+ */
+void msvmmaj_make_kernel(struct MajModel *model, struct MajData *data)
+{
+	if (model->kerneltype == K_LINEAR)
+		return;
+
+	long i, j;
+	long n = model->n;
+	double value;
+	double *x1, *x2;
+	double *K = Calloc(double, n*n*sizeof(double));
+
+	for (i=0; i<n; i++) {
+		for (j=i; j<n; j++) {
+			x1 = &data->Z[i*(data->m+1)+1];
+			x2 = &data->Z[j*(data->m+1)+1];
+			if (model->kerneltype == K_POLY)
+				value = msvmmaj_compute_poly(x1, x2, 
+						model->kernelparam, data->m);
+			else if (model->kerneltype == K_RBF)
+				value = msvmmaj_compute_rbf(x1, x2, 
+						model->kernelparam, data->m);
+			else if (model->kerneltype == K_SIGMOID)
+				value = msvmmaj_compute_rbf(x1, x2, 
+						model->kernelparam, data->m);
+			else {
+				fprintf(stderr, "Unknown kernel type in "
+						"msvmmaj_make_kernel\n");
+				exit(1);
+			}
+			matrix_set(K, n, i, j, value);
+			matrix_set(K, n, j, i, value);
+		}
+	}
+
+	// get cholesky if necessary.
+	if (model->use_cholesky == true) {
+		int status = dpotrf('L', n, K, n);
+		if (status != 0) {
+			fprintf(stderr, "Error (%i) computing Cholesky "
+					"decomposition of kernel matrix.\n",
+					status);
+			exit(0);
+		}
+		note("Got Cholesky.\n");
+	}
+
+	// copy kernel/cholesky to data
+	data->Z = realloc(data->Z, n*(n+1)*(sizeof(double)));
+	for (i=0; i<n; i++) {
+		for (j=0; j<n; j++)
+			matrix_set(data->Z, n+1, i, j+1, 
+					matrix_get(K, n, i, j));
+		matrix_set(data->Z, n+1, i, 0, 1.0);
+	}
+	data->m = n;
+	
+	// let data know what it's made of
+	data->kerneltype = model->kerneltype;
+	free(data->kernelparam);
+	switch (model->kerneltype) {
+		case K_LINEAR:
+			break;
+		case K_POLY:
+			data->kernelparam = Calloc(double, 3);
+			data->kernelparam[0] = model->kernelparam[0];
+			data->kernelparam[1] = model->kernelparam[1];
+			data->kernelparam[2] = model->kernelparam[2];
+			break;
+		case K_RBF:
+			data->kernelparam = Calloc(double, 1);
+			data->kernelparam[0] = model->kernelparam[0];
+			break;
+		case K_SIGMOID:
+			data->kernelparam = Calloc(double, 2);
+			data->kernelparam[0] = model->kernelparam[0];
+			data->kernelparam[1] = model->kernelparam[1];
+	}
+	data->use_cholesky = model->use_cholesky;
+	model->m = n;
+	free(K);
+}
+
+/**
+ * @brief Compute the RBF kernel between two vectors
+ * 
+ * @details
+ * The RBF kernel is computed between two vectors. This kernel is defined as
+ * @f[
+ * 	k(x_1, x_2) = \exp( -\gamma \| x_1 - x_2 \|^2 )
+ * @f]
+ * where @f$ \gamma @f$ is a kernel parameter specified.
+ *
+ * @param[in] 	x1 		first vector
+ * @param[in] 	x2 		second vector
+ * @param[in] 	kernelparam 	array of kernel parameters (gamma is first 
+ * 				element)
+ * @param[in] 	n 		length of the vectors x1 and x2
+ * @returns  			kernel evaluation
+ */
+double msvmmaj_compute_rbf(double *x1, double *x2, double *kernelparam, long n)
+{
+	long i;
+	double value = 0.0;
+	
+	for (i=0; i<n; i++) 
+		value += (x1[i] - x2[i]) * (x1[i] - x2[i]);
+	value *= -kernelparam[0];
+	return exp(value);
+}
+
+/**
+ * @brief Compute the polynomial kernel between two vectors
+ *
+ * @details
+ * The polynomial kernel is computed between two vectors. This kernel is 
+ * defined as
+ * @f[
+ * 	k(x_1, x_2) = ( \gamma \langle x_1, x_2 \rangle + c)^d
+ * @f]
+ * where @f$ \gamma @f$, @f$ c @f$ and @f$ d @f$ are kernel parameters.
+ *
+ * @param[in] 	x1 		first vector
+ * @param[in] 	x2 		second vector
+ * @param[in] 	kernelparam 	array of kernel parameters (gamma, c, d)
+ * @param[in] 	n 		length of the vectors x1 and x2
+ * @returns 			kernel evaluation
+ */
+double msvmmaj_compute_poly(double *x1, double *x2, double *kernelparam, long n)
+{
+	long i;
+	double value = 0.0;
+	for (i=0; i<n; i++)
+		value += x1[i]*x2[i];
+	value *= kernelparam[0];
+	value += kernelparam[1];
+	return pow(value, ((int) kernelparam[2]));
+}
+
+/**
+ * @brief Compute the sigmoid kernel between two vectors
+ * 
+ * @details
+ * The sigmoid kernel is computed between two vectors. This kernel is defined
+ * as
+ * @f[
+ * 	k(x_1, x_2) = \tanh( \gamma \langle x_1 , x_2 \rangle + c)
+ * @f]
+ * where @f$ \gamma @f$ and @f$ c @f$ are kernel parameters.
+ *
+ * @param[in] 	x1 		first vector
+ * @param[in] 	x2 		second vector
+ * @param[in] 	kernelparam 	array of kernel parameters (gamma, c)
+ * @param[in] 	n 		length of the vectors x1 and x2
+ * @returns 			kernel evaluation
+ */
+double msvmmaj_compute_sigmoid(double *x1, double *x2, double *kernelparam, long n)
+{
+	long i;
+	double value = 0.0;
+	for (i=0; i<n; i++)
+		value += x1[i]*x2[i];
+	value *= kernelparam[0];
+	value += kernelparam[1];
+	return tanh(value);
+}
diff --git a/src/msvmmaj_lapack.c b/src/msvmmaj_lapack.c
new file mode 100644
index 0000000..9ca8dab
--- /dev/null
+++ b/src/msvmmaj_lapack.c
@@ -0,0 +1,129 @@
+/**
+ * @file msvmmaj_lapack.c
+ * @author Gertjan van den Burg
+ * @date August 9, 2013
+ * @brief Utility functions for interacting with LAPACK
+ *
+ * @details
+ * Functions in this file are auxiliary functions which make it easier
+ * to use LAPACK functions from liblapack.
+ */
+
+#include "msvmmaj_lapack.h"
+
+/**
+ * @brief Solve AX = B where A is symmetric positive definite.
+ *
+ * @details
+ * Solve a linear system of equations AX = B where A is symmetric positive 
+ * definite. This function uses the externel LAPACK routine dposv.
+ *
+ * @param[in] 		UPLO 	which triangle of A is stored
+ * @param[in] 		N 	order of A
+ * @param[in] 		NRHS 	number of columns of B
+ * @param[in,out] 	A 	double precision array of size (LDA, N). On 
+ * 				exit contains the upper or lower factor of the
+ * 				Cholesky factorization of A.
+ * @param[in] 		LDA 	leading dimension of A
+ * @param[in,out] 	B 	double precision array of size (LDB, NRHS). On
+ * 				exit contains the N-by-NRHS solution matrix X.
+ * @param[in] 		LDB 	the leading dimension of B
+ * @returns 			info parameter which contains the status of the
+ * 				computation:
+ * 					- =0: 	success
+ * 					- <0: 	if -i, the i-th argument had
+ * 						an illegal value
+ * 					- >0: 	if i, the leading minor of A
+ * 						was not positive definite
+ *
+ * See the LAPACK documentation at: 
+ * http://www.netlib.org/lapack/explore-html/dc/de9/group__double_p_osolve.html
+ */
+int dposv(char UPLO, int N, int NRHS, double *A, int LDA, double *B,
+		int LDB)
+{
+	extern void dposv_(char *UPLO, int *Np, int *NRHSp, double *A,
+			int *LDAp, double *B, int *LDBp, int *INFOp);
+	int INFO;
+	dposv_(&UPLO, &N, &NRHS, A, &LDA, B, &LDB, &INFO);
+	return INFO;
+}
+
+/**
+ * @brief Solve a system of equations AX = B where A is symmetric.
+ *
+ * @details
+ * Solve a linear system of equations AX = B where A is symmetric. This 
+ * function uses the external LAPACK routine dsysv.
+ *
+ * @param[in] 		UPLO 	which triangle of A is stored
+ * @param[in] 		N 	order of A
+ * @param[in] 		NRHS 	number of columns of B
+ * @param[in,out] 	A 	double precision array of size (LDA, N). On 
+ * 				exit contains the block diagonal matrix D and
+ * 				the multipliers used to obtain the factor U or
+ * 				L from the factorization A = U*D*U**T or
+ * 				A = L*D*L**T.
+ * @param[in] 		LDA 	leading dimension of A
+ * @param[in] 		IPIV 	integer array containing the details of D
+ * @param[in,out] 	B 	double precision array of size (LDB, NRHS). On
+ * 				exit contains the N-by-NRHS matrix X
+ * @param[in] 		LDB 	leading dimension of B
+ * @param[out] 		WORK 	double precision array of size max(1,LWORK). On
+ * 				exit, WORK(1) contains the optimal LWORK
+ * @param[in] 		LWORK 	the length of WORK, can be used for determining
+ * 				the optimal blocksize for dsystrf.
+ * @returns 			info parameter which contains the status of the
+ * 				computation:
+ * 					- =0: 	success
+ * 					- <0: 	if -i, the i-th argument had an
+ * 						illegal value
+ * 					- >0: 	if i, D(i, i) is exactly zero,
+ * 						no solution can be computed.
+ *
+ * See the LAPACK documentation at:
+ * http://www.netlib.org/lapack/explore-html/d6/d0e/group__double_s_ysolve.html
+ */
+int dsysv(char UPLO, int N, int NRHS, double *A, int LDA, int *IPIV,
+		double *B, int LDB, double *WORK, int LWORK)
+{
+	extern void dsysv_(char *UPLO, int *Np, int *NRHSp, double *A,
+			int *LDAp, int *IPIV, double *B, int *LDBp,
+			double *WORK, int *LWORK, int *INFOp);
+	int INFO;
+	dsysv_(&UPLO, &N, &NRHS, A, &LDA, IPIV, B, &LDB, WORK, &LWORK, &INFO);
+	return INFO;
+}
+
+/**
+ * @brief Compute the Cholesky factorization of a real symmetric positive 
+ *  	  definite matrix.
+ *
+ * @details
+ * This function uses the external LAPACK routine dpotrf.
+ *
+ * @param[in] 		UPLO 	which triangle of A is stored
+ * @param[in] 		N 	order of A
+ * @param[in,out] 	A 	double precision array of size (LDA, N). On
+ * 				exit contains the factor U or L of the Cholesky
+ * 				factorization
+ * @param[in] 		LDA 	leading dimension of A
+ * @returns 			info parameter which contains the status of the
+ * 				computation:
+ * 					- =0: 	success
+ * 					- <0: 	if -i, the i-th argument had an
+ * 						illegal value
+ * 					- >0: 	if i, the leading minor of
+ * 						order i is not positive
+ * 						definite
+ *
+ * See the LAPACK documentation at:
+ * http://www.netlib.org/lapack/explore-html/d0/d8a/dpotrf_8f.html
+ */
+int dpotrf(char UPLO, int N, double *A, int LDA)
+{
+	extern void dpotrf_(char *UPLO, int *N, double *A, int *LDA, int *INFOp);
+	int INFO;
+	dpotrf_(&UPLO, &N, A, &LDA, &INFO);
+	return INFO;
+}
diff --git a/src/msvmmaj_matrix.c b/src/msvmmaj_matrix.c
new file mode 100644
index 0000000..ffa0c21
--- /dev/null
+++ b/src/msvmmaj_matrix.c
@@ -0,0 +1,153 @@
+/**
+ * @file msvmmaj_matrix.c
+ * @author Gertjan van den Burg
+ * @date August, 2013
+ * @brief Functions facilitating matrix access
+ *
+ * @details
+ * The functions contained in this file are used when
+ * accessing or writing to matrices. Seperate functions
+ * exist of adding and multiplying existing matrix
+ * elements, to ensure this is done in place.
+ *
+ */
+
+#include "msvmmaj_matrix.h"
+#include "util.h"
+
+/**
+ * @brief Set element of matrix 
+ *
+ * @details
+ * Row-Major order is used to set a matrix element. Since matrices
+ * of type double are most common in MSVMMaj, this function only
+ * deals with that type.
+ *
+ * @param[in] 	M 	matrix to set element of
+ * @param[in] 	cols 	number of columns of M
+ * @param[in] 	i 	row index of element to write to
+ * @param[in] 	j 	column index of element to write to
+ * @param[out] 	val 	value to write to specified element of M
+ */
+void matrix_set(double *M, long cols, long i, long j, double val)
+{
+	M[i*cols+j] = val;
+}
+
+/**
+ * @brief Retrieve value from matrix
+ * 
+ * @details
+ * Return a value from a matrix using row-major order.
+ *
+ * @param[in] 	M 	matrix to retrieve value from
+ * @param[in] 	cols 	number of columns of M
+ * @param[in] 	i 	row index (starting from 0)
+ * @param[in] 	j 	column index (starting from 0)
+ * @returns 		matrix element at (i, j)	
+ */
+double matrix_get(double *M, long cols, long i, long j)
+{
+	return M[i*cols+j];
+}
+
+/**
+ * @brief Add value to matrix element
+ *
+ * @details
+ * This function is added to efficiently add values to matrix
+ * elements, without having to use get and set methods.
+ *
+ * @param[in] 	M 	matrix
+ * @param[in] 	cols 	number of columns of M
+ * @param[in] 	i  	row index (starting from 0)
+ * @param[in] 	j 	column index (starting from 0)
+ * @param[in] 	val 	value to add to matrix element (i, j)
+ */
+void matrix_add(double *M, long cols, long i, long j, double val)
+{
+	M[i*cols+j] += val;
+}
+
+/**
+ * @brief Multiply matrix element by value
+ * 
+ * @details
+ * This function is added to efficiently multiply a matrix element
+ * by a certain value, without having to use get and set methods.
+ *
+ * @param[in] 	M 	matrix
+ * @param[in] 	cols 	number of columns of M
+ * @param[in] 	i 	row index (starting from 0)
+ * @param[in] 	j 	column index (starting from 0)
+ * @param[in] 	val 	value to multiply matrix element (i, j) with
+ */
+void matrix_mul(double *M, long cols, long i, long j, double val)
+{
+	M[i*cols+j] *= val;
+}
+
+/**
+ * @brief Set element of 3D matrix
+ *
+ * @details
+ * Set an element of a 3D matrix using row-major order. 
+ *
+ * @param[in] 	M 	matrix
+ * @param[in] 	N2 	second dimension of M
+ * @param[in] 	N3 	third dimension of M
+ * @param[in] 	i 	index along first dimension
+ * @param[in] 	j 	index along second dimension
+ * @param[in] 	k 	index along third dimension
+ * @param[in] 	val 	value to set element (i, j, k) to
+ *
+ * See:
+ * http://en.wikipedia.org/wiki/Row-major_order
+ */
+void matrix3_set(double *M, long N2, long N3, long i, long j, 
+		long k, double val)
+{
+	M[k+N3*(j+N2*i)] = val;
+}
+
+/**
+ * @brief Get element of 3D matrix
+ *
+ * @details
+ * Retrieve an element from a 3D matrix.
+ *
+ * @param[in] 	M 	matrix
+ * @param[in] 	N2 	second dimension of M
+ * @param[in] 	N3 	third dimension of M
+ * @param[in] 	i 	index along first dimension
+ * @param[in] 	j 	index along second dimension
+ * @param[in] 	k 	index along third dimension
+ * @returns 		value at the (i, j, k) element of M
+ */
+double matrix3_get(double *M, long N2, long N3, long i, long j,
+		long k)
+{
+	return M[k+N3*(j+N2*i)];
+}
+
+/**
+ * @brief print a matrix
+ * 
+ * @details
+ * Debug function to print a matrix
+ *
+ * @param[in] 	M 	matrix
+ * @param[in] 	rows 	number of rows of M
+ * @param[in] 	cols 	number of columns of M
+ */
+void print_matrix(double *M, long rows, long cols)
+{
+	long i, j;
+
+	for (i=0; i<rows; i++) {
+		for (j=0; j<cols; j++)
+			note("%8.8f ", matrix_get(M, cols, i, j));
+		note("\n");
+	}
+	note("\n");
+}
diff --git a/src/msvmmaj_pred.c b/src/msvmmaj_pred.c
index 5f1b1ae..98b6e0a 100644
--- a/src/msvmmaj_pred.c
+++ b/src/msvmmaj_pred.c
@@ -1,31 +1,36 @@
 /**
  * @file msvmmaj_pred.c
- * @author Gertjan van den Burg (burg@ese.eur.nl)
+ * @author Gertjan van den Burg
  * @date August 9, 2013
  * @brief Main functions for predicting class labels..
  *
+ * @details
+ * This file contains functions for predicting the class labels of instances
+ * and a function for calculating the predictive performance (hitrate) of 
+ * a prediction given true class labels.
+ *
  */
 
 #include <cblas.h>
 
 #include "libMSVMMaj.h"
-#include "MSVMMaj.h"
-#include "matrix.h"
+#include "msvmmaj.h"
+#include "msvmmaj_matrix.h"
 #include "msvmmaj_pred.h"
 
 /**
- * @name predict_labels
  * @brief Predict class labels of data given and output in predy
  *
+ * @details
  * The labels are predicted by mapping each instance in data to the 
  * simplex space using the matrix V in the given model. Next, for each
  * instance the nearest simplex vertex is determined using an Euclidean
  * norm. The nearest simplex vertex determines the predicted class label,
- * which is recorded in predy
+ * which is recorded in predy.
  *
- * @param [in] 		data 		data to predict labels for
- * @param [in] 		model 		model with optimized V
- * @param [out] 	predy 		pre-allocated vector to record predictions in
+ * @param[in] 	data 		MajData to predict labels for
+ * @param[in] 	model 		MajModel with optimized V
+ * @param[out] 	predy 		pre-allocated vector to record predictions in
  */
 void msvmmaj_predict_labels(struct MajData *data, struct MajModel *model, long *predy)
 {
@@ -84,15 +89,15 @@ void msvmmaj_predict_labels(struct MajData *data, struct MajModel *model, long *
 }
 
 /**
- * @name msvmmaj_prediction_perf
  * @brief Calculate the predictive performance (percentage correct)
  *
+ * @details
  * The predictive performance is calculated by simply counting the number
  * of correctly classified samples and dividing by the total number of 
  * samples, multiplying by 100.
  *
- * @param [in] 		data 		the dataset with known labels
- * @param [in] 		predy 		the predicted class labels
+ * @param[in] 	data 	the MajData dataset with known labels
+ * @param[in] 	predy 	the predicted class labels
  *
  * @returns percentage correctly classified.
  */
diff --git a/src/msvmmaj_train.c b/src/msvmmaj_train.c
index 272d86a..97ee6a1 100644
--- a/src/msvmmaj_train.c
+++ b/src/msvmmaj_train.c
@@ -1,6 +1,6 @@
 /**
  * @file msvmmaj_train.c
- * @author Gertjan van den Burg (burg@ese.eur.nl)
+ * @author Gertjan van den Burg
  * @date August 9, 2013
  * @brief Main functions for training the MSVMMaj solution.
  * 
@@ -13,25 +13,34 @@
 #include <math.h>
 #include <cblas.h>
 
-#include "msvmmaj_train.h"
-#include "MSVMMaj.h"
 #include "libMSVMMaj.h"
-#include "mylapack.h"
-#include "matrix.h"
+#include "msvmmaj.h"
+#include "msvmmaj_lapack.h"
+#include "msvmmaj_matrix.h"
+#include "msvmmaj_train.h"
 #include "util.h"
 
+/**
+ * Maximum number of iterations of the algorithm.
+ */
 #define MAX_ITER 1000000
 
 /**
- * @name msvmmaj_optimize
  * @brief The main training loop for MSVMMaj
  *
- * The msvmmaj_optimize() function is the main training function. This function
+ * @details
+ * This function is the main training function. This function
  * handles the optimization of the model with the given model parameters, with
- * the data given. On return the matrix model->V contains the optimal weight matrix.
+ * the data given. On return the matrix MajModel::V contains the optimal 
+ * weight matrix.
+ *
+ * In this function, step doubling is used in the majorization algorithm after
+ * a burn-in of 50 iterations. If the training is finished, MajModel::t and 
+ * MajModel::W are extracted from MajModel::V.
  *
- * @param [in,out] 	model 	the model to be trained. Contains optimal V on exit.
- * @param [in] 		data 	the data to train the model with.
+ * @param[in,out] 	model 	the MajModel to be trained. Contains optimal 
+ * 				V on exit.
+ * @param[in] 		data 	the MajData to train the model with.
  */
 void msvmmaj_optimize(struct MajModel *model, struct MajData *data)
 {
@@ -49,7 +58,7 @@ void msvmmaj_optimize(struct MajModel *model, struct MajData *data)
 	double *ZAZVT = Calloc(double, (m+1)*(K-1));
 
 	note("Starting main loop.\n");
-	note("MajDataset:\n");
+	note("Dataset:\n");
 	note("\tn = %i\n", n);
 	note("\tm = %i\n", m);
 	note("\tK = %i\n", K);
@@ -78,8 +87,8 @@ void msvmmaj_optimize(struct MajModel *model, struct MajData *data)
 		L = msvmmaj_get_loss(model, data, ZV);
 
 		if (it%50 == 0)
-			note("iter = %li, L = %15.16f, Lbar = %15.16f, reldiff = %15.16f\n",
-					it, L, Lbar, (Lbar - L)/L);
+			note("iter = %li, L = %15.16f, Lbar = %15.16f, "
+			     "reldiff = %15.16f\n", it, L, Lbar, (Lbar - L)/L);
 		it++;
 	}
 
@@ -91,7 +100,8 @@ void msvmmaj_optimize(struct MajModel *model, struct MajData *data)
 		model->t[i] = matrix_get(model->V, K-1, 0, i);
 	for (i=1; i<m+1; i++)
 		for (j=0; j<K-1; j++)
-			matrix_set(model->W, K-1, i-1, j, matrix_get(model->V, K-1, i, j));
+			matrix_set(model->W, K-1, i-1, j, 
+					matrix_get(model->V, K-1, i, j));
 	free(B);
 	free(ZV);
 	free(ZAZ);
@@ -100,19 +110,22 @@ void msvmmaj_optimize(struct MajModel *model, struct MajData *data)
 }
 
 /**
- * @name msvmmaj_get_loss
- * @brief calculate the current value of the loss function 
+ * @brief Calculate the current value of the loss function
  * 
- * The current loss value is calculated based on the matrix V in the given
- * model.
+ * @details
+ * The current loss function value is calculated based on the matrix V in the 
+ * given model. Note that the matrix ZV is passed explicitly to avoid having
+ * to reallocate memory at every step.
  *
- * @param [in] 		model 	model structure which holds the current estimate V
- * @param [in]  	data 	data structure
- * @param [in,out] 	ZV 	pre-allocated matrix ZV which is updated on output
- * 
- * @return 	the current value of the loss function
+ * @param[in] 		model 	MajModel structure which holds the current 
+ * 				estimate V
+ * @param[in]  		data 	MajData structure
+ * @param[in,out] 	ZV 	pre-allocated matrix ZV which is updated on 
+ * 				output
+ * @returns 			the current value of the loss function
  */
-double msvmmaj_get_loss(struct MajModel *model, struct MajData *data, double *ZV)
+double msvmmaj_get_loss(struct MajModel *model, struct MajData *data, 
+		double *ZV)
 {
 	long i, j;
 	long n = data->n;
@@ -151,10 +164,52 @@ double msvmmaj_get_loss(struct MajModel *model, struct MajData *data, double *ZV
 }
 
 /**
- * @name msvmmaj_get_update
- * @brief perform a single step of the majorization algorithm to update V
+ * @brief Perform a single step of the majorization algorithm to update V
+ *
+ * @details
+ * This function contains the main update calculations of the algorithm. These
+ * calculations are necessary to find a new update V. The calculations exist of
+ * recalculating the majorization coefficients for all instances and all
+ * classes, and solving a linear system to find V.
+ *
+ * Because the function msvmmaj_get_update() is always called after a call to
+ * msvmmaj_get_loss() with the same MajModel::V, it is unnecessary to calculate
+ * the updated errors MajModel::Q and MajModel::H here too. This saves on 
+ * computation time.
  *
- * details
+ * In calculating the majorization coefficients we calculate the elements of a 
+ * diagonal matrix A with elements
+ * @f[
+ * 	A_{i, i} = \frac{1}{n} \rho_i \sum_{j \neq k} \left[ 
+ * 		\varepsilon_i a_{ijk}^{(p)} + (1 - \varepsilon_i) \omega_i
+ * 		a_{ijk}^{(p)} \right],
+ * @f]
+ * where @f$ k = y_i @f$.
+ * Since this matrix is only used to calculate the matrix @f$ Z' A Z @f$, it is
+ * efficient to update a matrix ZAZ through consecutive rank 1 updates with 
+ * a single element of A and the corresponding row of Z. The BLAS function
+ * dsyr is used for this.
+ *
+ * The B matrix is has rows
+ * @f[
+ * 	\boldsymbol{\beta}_i' = \frac{1}{n} \rho_i \sum_{j \neq k} \left[
+ * 		\varepsilon_i \left( b_{ijk}^{(1)} - a_{ijk}^{(1)}
+ * 			\overline{q}_i^{(kj)} \right) + (1 - \varepsilon_i)
+ * 		\omega_i \left( b_{ijk}^{(p)} - a_{ijk}^{(p)}
+ * 			\overline{q}_i^{(kj)} \right) \right]
+ * 		\boldsymbol{\delta}_{kj}'
+ * @f]
+ * This is also split into two cases, one for which @f$ \varepsilon_i = 1 @f$,
+ * and one for when it is 0. The 3D simplex difference matrix is used here, in
+ * the form of the @f$ \boldsymbol{\delta}_{kj}' @f$.
+ *
+ * Finally, the following system is solved
+ * @f[
+ * 	(\textbf{Z}'\textbf{AZ} + \lambda \textbf{J})\textbf{V} =
+ * 		(\textbf{Z}'\textbf{AZ}\overline{\textbf{V}} + \textbf{Z}'
+ * 		\textbf{B})
+ * @f]
+ * solving this system is done through dposv().
  *
  * @param [in,out] 	model 	model to be updated
  * @param [in] 		data 	data used in model
@@ -166,9 +221,6 @@ double msvmmaj_get_loss(struct MajModel *model, struct MajData *data, double *ZV
 void msvmmaj_get_update(struct MajModel *model, struct MajData *data, double *B,
 		double *ZAZ, double *ZAZV, double *ZAZVT)
 {
-	// Because msvmmaj_update is always called after a call to 
-	// msvmmaj_get_loss() with the latest V, it is unnecessary to recalculate
-	// the matrix ZV, the errors Q, or the Huber errors H. Awesome!
 	int status, class;
 	long i, j, k;
 	double Avalue, Bvalue;
@@ -182,11 +234,14 @@ void msvmmaj_get_update(struct MajModel *model, struct MajData *data, double *B,
 	double p = model->p;
 	double *rho = model->rho;
 
+	// constants which are used often throughout
 	const double a2g2 = 0.25*p*(2.0*p - 1.0)*pow((kappa+1.0)/2.0,p-2.0);
 	const double in = 1.0/((double) n);
 
+	// clear matrices
 	Memset(B, double, n*(K-1));
 	Memset(ZAZ, double, (m+1)*(m+1));
+
 	b = 0;
 	for (i=0; i<n; i++) {
 		value = 0;
@@ -215,7 +270,8 @@ void msvmmaj_get_update(struct MajModel *model, struct MajData *data, double *B,
 					b = 0;
 				}
 				for (k=0; k<K-1; k++) {
-					Bvalue = in*rho[i]*b*matrix3_get(model->UU, K-1, K, i, k, j);
+					Bvalue = in*rho[i]*b*matrix3_get(
+						model->UU, K-1, K, i, k, j);
 					matrix_add(B, K-1, i, k, Bvalue);
 				}
 				Avalue += a*matrix_get(model->R, K, i, j);
@@ -227,13 +283,27 @@ void msvmmaj_get_update(struct MajModel *model, struct MajData *data, double *B,
 					if (q <= -kappa) {
 						b = 0.5 - kappa/2.0 - q;
 					} else if ( q <= 1.0) {
-						b = pow(1.0 - q, 3.0)/(2.0*pow(kappa + 1.0, 2.0));
+						b = pow(1.0 - q, 3.0)/(
+							2.0*pow(kappa + 1.0, 
+								2.0));
 					} else {
 						b = 0;
 					}
 					for (k=0; k<K-1; k++) {
-						Bvalue = in*rho[i]*omega*b*matrix3_get(model->UU, K-1, K, i, k, j);
-						matrix_add(B, K-1, i, k, Bvalue);
+						Bvalue = in*rho[i]*omega*b*
+							matrix3_get(
+								model->UU,
+								K-1,
+							       	K,
+								i,
+								k,
+								j);
+						matrix_add(
+								B,
+								K-1,
+								i,
+								k,
+								Bvalue);
 					}
 				}
 				Avalue = 1.5*(K - 1.0);
@@ -241,23 +311,51 @@ void msvmmaj_get_update(struct MajModel *model, struct MajData *data, double *B,
 				for (j=0; j<K; j++) {
 					q = matrix_get(model->Q, K, i, j);
 					if (q <= (p + kappa - 1.0)/(p - 2.0)) {
-						a = 0.25*pow(p, 2.0)*pow(0.5 - kappa/2.0 - q, p - 2.0);
+						a = 0.25*pow(p, 2.0)*pow(
+							0.5 - kappa/2.0 - q,
+						       		p - 2.0);
 					} else if (q <= 1.0) {
 						a = a2g2;
 					} else {
-						a = 0.25*pow(p, 2.0)*pow((p/(p - 2.0))*(0.5 - kappa/2.0 - q), p - 2.0);
-						b = a*(2.0*q + kappa - 1.0)/(p - 2.0) + 0.5*p*pow((p/(p - 2.0))*(0.5 - kappa/2.0 - q), p - 1.0);
+						a = 0.25*pow(p, 2.0)*pow(
+							(p/(p - 2.0))*
+							(0.5 - kappa/2.0 - q), 
+							p - 2.0);
+						b = a*(2.0*q + kappa - 1.0)/
+							(p - 2.0) + 
+							0.5*p*pow(
+								p/(p - 2.0)*
+								(0.5 - kappa/
+								 2.0 - q), 
+								p - 1.0);
 					}
 					if (q <= -kappa) {
-						b = 0.5*p*pow(0.5 - kappa/2.0 - q, p - 1.0);
+						b = 0.5*p*pow(
+							0.5 - kappa/2.0 - q, 
+							p - 1.0);
 					} else if ( q <= 1.0) {
-						b = p*pow(1.0 - q, 2.0*p - 1.0)/pow(2*kappa+2.0, p);
+						b = p*pow(1.0 - q, 
+								2.0*p - 1.0)/
+							pow(2*kappa+2.0, p);
 					}
 					for (k=0; k<K-1; k++) {
-						Bvalue = in*rho[i]*omega*b*matrix3_get(model->UU, K-1, K, i, k, j);
-						matrix_add(B, K-1, i, k, Bvalue);
+						Bvalue = in*rho[i]*omega*b*
+							matrix3_get(
+								model->UU,
+								K-1,
+								K,
+								i,
+								k,
+								j);
+						matrix_add(
+								B,
+								K-1,
+								i,
+								k,
+								Bvalue);
 					}
-					Avalue += a*matrix_get(model->R, K, i, j);
+					Avalue += a*matrix_get(model->R,
+						       	K, i, j);
 				}
 			}
 			Avalue *= omega;
@@ -352,7 +450,8 @@ void msvmmaj_get_update(struct MajModel *model, struct MajData *data, double *B,
 		// This step should not be necessary, as the matrix
 		// ZAZ is positive semi-definite by definition. It 
 		// is included for safety.
-		fprintf(stderr, "Received nonzero status from dposv: %i\n", status);
+		fprintf(stderr, "Received nonzero status from dposv: %i\n", 
+				status);
 		int *IPIV = malloc((m+1)*sizeof(int));
 		double *WORK = malloc(1*sizeof(double));
 		status = dsysv(
@@ -379,7 +478,8 @@ void msvmmaj_get_update(struct MajModel *model, struct MajData *data, double *B,
 				WORK,
 				sizeof(WORK)/sizeof(double));
 		if (status != 0)
-			fprintf(stderr, "Received nonzero status from dsysv: %i\n", status);
+			fprintf(stderr, "Received nonzero status from "
+					"dsysv: %i\n", status);
 	}
 
 	// Return to Row-major order. The matrix ZAZVT contains the 
@@ -403,8 +503,18 @@ void msvmmaj_get_update(struct MajModel *model, struct MajData *data, double *B,
 	
 	for (i=0; i<m+1; i++) {
 		for (j=0; j<K-1; j++) {
-			matrix_set(model->Vbar, K-1, i, j, matrix_get(model->V, K-1, i, j));
-			matrix_set(model->V, K-1, i, j, matrix_get(ZAZV, K-1, i, j));
+			matrix_set(
+					model->Vbar,
+					K-1,
+					i,
+					j,
+					matrix_get(model->V, K-1, i, j));
+			matrix_set(
+					model->V,
+					K-1,
+					i,
+					j,
+					matrix_get(ZAZV, K-1, i, j));
 		}
 	}
 }
diff --git a/src/msvmmaj_train_dataset.c b/src/msvmmaj_train_dataset.c
index 2da8bee..4f5f4d9 100644
--- a/src/msvmmaj_train_dataset.c
+++ b/src/msvmmaj_train_dataset.c
@@ -1,22 +1,53 @@
+/**
+ * @file msvmmaj_train_dataset.c
+ * @author Gertjan van den Burg
+ * @date January, 2014
+ * @brief Functions for finding the optimal parameters for the dataset
+ *
+ * @details
+ * The MSVMMaj algorithm takes a number of parameters. The functions in
+ * this file are used to find the optimal parameters. 
+ */
+
 #include <math.h>
 #include <time.h>
 
 #include "crossval.h"
 #include "libMSVMMaj.h"
-#include "matrix.h"
+#include "msvmmaj.h"
+#include "msvmmaj_init.h"
+#include "msvmmaj_matrix.h"
 #include "msvmmaj_train.h"
 #include "msvmmaj_train_dataset.h"
 #include "msvmmaj_pred.h"
-#include "MSVMMaj.h"
 #include "util.h"
 #include "timer.h"
 
 extern FILE *MSVMMAJ_OUTPUT_FILE;
 
+/**
+ * @brief Initialize a Queue from a Training instance
+ *
+ * @details
+ * A Training instance describes the grid to search over. This funtion
+ * creates all tasks that need to be performed and adds these to 
+ * a Queue. Each task contains a pointer to the train and test datasets
+ * which are supplied. Note that the tasks are created in a specific order of
+ * the parameters, to ensure that the MajModel::V of a previous parameter
+ * set provides the best possible initial estimate of MajModel::V for the next
+ * parameter set.
+ *
+ * @param[in] 	training 	Training struct describing the grid search
+ * @param[in] 	queue 		pointer to a Queue that will be used to 
+ * 				add the tasks to
+ * @param[in] 	train_data 	MajData of the training set
+ * @param[in] 	test_data 	MajData of the test set
+ *
+ */
 void make_queue(struct Training *training, struct Queue *queue, 
 		struct MajData *train_data, struct MajData *test_data)
 {
-	long i, j, k, l, m;
+	long i, j, k;
 	long N, cnt = 0;
 	struct Task *task;
 	queue->i = 0;
@@ -26,30 +57,122 @@ void make_queue(struct Training *training, struct Queue *queue,
 	N *= training->Nk;
 	N *= training->Ne;
 	N *= training->Nw;
+	// these parameters are not necessarily non-zero
+	N *= training->Ng > 0 ? training->Ng : 1;
+	N *= training->Nc > 0 ? training->Nc : 1;
+	N *= training->Nd > 0 ? training->Nd : 1;
 
 	queue->tasks = Malloc(struct Task *, N);
 	queue->N = N;
 
-	for (i=0; i<training->Ne; i++)
+	// initialize all tasks
+	for (i=0; i<N; i++) {
+		task = Malloc(struct Task, 1);
+		task->ID = i;
+		task->train_data = train_data;
+		task->test_data = test_data;
+		task->folds = training->folds;
+		task->kerneltype = training->kerneltype;
+		task->kernel_param = Calloc(double, training->Ng + 
+				training->Nc + training->Nd);
+		queue->tasks[i] = task;
+	}
+
+	// These loops mimick a large nested for loop. The advantage is that
+	// Nd, Nc and Ng which are on the outside of the nested for loop can
+	// now be zero, without large modification (see below). Whether this
+	// is indeed better than the nested for loop has not been tested.
+	cnt = 1;
+	i = 0;
+	while (i < N )
+		for (j=0; j<training->Np; j++) 
+			for (k=0; k<cnt; k++) {
+				queue->tasks[i]->p = training->ps[j];
+				i++;
+			}
+
+	cnt *= training->Np;
+	i = 0;
+	while (i < N )
+		for (j=0; j<training->Nl; j++) 
+			for (k=0; k<cnt; k++) {
+				queue->tasks[i]->lambda = 
+					training->lambdas[j];
+				i++;
+			}
+
+	cnt *= training->Nl;
+	i = 0;
+	while (i < N )
+		for (j=0; j<training->Nk; j++)
+			for (k=0; k<cnt; k++) {
+				queue->tasks[i]->kappa = training->kappas[j];
+				i++;
+			}
+
+	cnt *= training->Nk;
+	i = 0;
+	while (i < N )
 		for (j=0; j<training->Nw; j++)
-			for (k=0; k<training->Nk; k++)
-				for (l=0; l<training->Nl; l++)
-					for (m=0; m<training->Np; m++) {
-						task = Malloc(struct Task, 1);
-						task->epsilon = training->epsilons[i];
-						task->weight_idx = training->weight_idxs[j];
-						task->kappa = training->kappas[k];
-						task->lambda = training->lambdas[l];
-						task->p = training->ps[m];
-						task->train_data = train_data;
-						task->test_data = test_data;
-						task->folds = training->folds;
-						task->ID = cnt;
-						queue->tasks[cnt] = task;
-						cnt++;
-					}
+			for (k=0; k<cnt; k++) {
+				queue->tasks[i]->weight_idx = 
+					training->weight_idxs[j];
+				i++;
+			}
+
+	cnt *= training->Nw;
+	i = 0;
+	while (i < N )
+		for (j=0; j<training->Ne; j++)
+			for (k=0; k<cnt; k++) {
+				queue->tasks[i]->epsilon = 
+					training->epsilons[j];
+				i++;
+			}
+
+	cnt *= training->Ne;
+	i = 0;
+	while (i < N && training->Ng > 0)
+		for (j=0; j<training->Ng; j++)
+			for (k=0; k<cnt; k++) {
+				queue->tasks[i]->kernel_param[0] = 
+					training->gammas[j];
+				i++;
+			}
+
+	cnt *= training->Ng > 0 ? training->Ng : 1;
+	i = 0;
+	while (i < N && training->Nc > 0)
+		for (j=0; j<training->Nc; j++)
+			for (k=0; k<cnt; k++) {
+				queue->tasks[i]->kernel_param[1] = 
+					training->coefs[j];
+				i++;
+			}
+
+	cnt *= training->Nc > 0 ? training->Ng : 1;
+	i = 0;
+	while (i < N && training->Nd > 0)
+		for (j=0; j<training->Nd; j++)
+			for (k=0; k<cnt; k++) {
+				queue->tasks[i]->kernel_param[2] = 
+					training->degrees[j];
+				i++;
+			}
 }
 
+/**
+ * @brief Get new Task from Queue
+ *
+ * @details
+ * Return a pointer to the next Task in the Queue. If no Task instances are
+ * left, NULL is returned. The internal counter Queue::i is used for finding
+ * the next Task.
+ *
+ * @param[in] 	q 	Queue instance
+ * @returns 		pointer to next Task
+ *
+ */
 struct Task *get_next_task(struct Queue *q)
 {
 	long i = q->i;
@@ -60,6 +183,19 @@ struct Task *get_next_task(struct Queue *q)
 	return NULL;
 }
 
+/**
+ * @brief Comparison function for Tasks based on performance
+ *
+ * @details
+ * To be able to sort Task structures on the performance of their specific
+ * set of parameters, this comparison function is implemented. Task structs
+ * are sorted with highest performance first.
+ *
+ * @param[in] 	elem1 	Task 1
+ * @param[in] 	elem2 	Task 2
+ * @returns 		result of inequality of Task 1 performance over
+ * 			Task 2 performance
+ */
 int tasksort(const void *elem1, const void *elem2)
 {
 	const struct Task *t1 = (*(struct Task **) elem1);
@@ -67,6 +203,16 @@ int tasksort(const void *elem1, const void *elem2)
 	return (t1->performance > t2->performance);
 }
 
+/**
+ * @brief Comparison function for doubles
+ *
+ * @details
+ * Similar to tasksort() only now for two doubles.
+ *
+ * @param[in] 	elem1 	number 1
+ * @param[in] 	elem2 	number 2
+ * @returns 		comparison of number 1 larger than number 2
+ */
 int doublesort(const void *elem1, const void *elem2)
 {
 	const double t1 = (*(double *) elem1);
@@ -74,7 +220,20 @@ int doublesort(const void *elem1, const void *elem2)
 	return t1 > t2;
 }
 
-
+/**
+ * @brief Calculate the percentile of an array of doubles
+ *
+ * @details
+ * The percentile of performance is used to find the top performing
+ * configurations. Since no standard definition of the percentile exists, we
+ * use the method used in MATLAB and Octave. Since calculating the percentile
+ * requires a sorted list of the values, a local copy is made first.
+ *
+ * @param[in] 	values 	array of doubles
+ * @param[in] 	N 	length of the array
+ * @param[in] 	p 	percentile to calculate ( 0 <= p <= 1.0 ).
+ * @returns 		the p-th percentile of the values
+ */
 double prctile(double *values, long N, double p)
 {
 	long i;
@@ -94,16 +253,50 @@ double prctile(double *values, long N, double p)
 	return boundary;
 }
 
+/**
+ * @brief Run repeats of the Task structs in Queue to find the best
+ * configuration
+ *
+ * @details
+ * The best performing tasks in the supplied Queue are found by taking those
+ * Task structs that have a performance greater or equal to the 95% percentile
+ * of the performance of all tasks. These tasks are then gathered in a new
+ * Queue. For each of the tasks in this new Queue the cross validation run is 
+ * repeated a number of times. 
+ *
+ * For each of the Task configurations that are repeated the mean performance,
+ * standard deviation of the performance and the mean computation time are
+ * reported. 
+ *
+ * Finally, the overall best tasks are written to the specified output. These
+ * tasks are selected to have both the highest mean performance, as well as the
+ * smallest standard deviation in their performance. This is done as follows.
+ * First the 99th percentile of task performance and the 1st percentile of
+ * standard deviation is calculated. If a task exists for which the mean
+ * performance of the repeats and the standard deviation equals these values
+ * respectively, this task is found to be the best and is written to the
+ * output. If no such task exists, the 98th percentile of performance and the
+ * 2nd percentile of standard deviation is considered. This is repeated until
+ * an interval is found which contains tasks. If one or more tasks are found,
+ * this loop stops.
+ *
+ * @param[in] 	q 		Queue of Task structs which have already been 
+ * 				run and have a Task::performance value
+ * @param[in] 	repeats 	Number of times to repeat the best
+ * 				configurations for consistency
+ * @param[in] 	traintype 	type of training to do (CV or TT)
+ *
+ */
 void consistency_repeats(struct Queue *q, long repeats, TrainType traintype)
 {
 	long i, r, N;
 	double p, pi, pr, boundary, time, *std, *mean, *perf;
 	struct Queue *nq = Malloc(struct Queue, 1);
-	struct MajModel *model = Malloc(struct MajModel, 1);
+	struct MajModel *model = msvmmaj_init_model();
 	struct Task *task = Malloc(struct Task, 1);
 	clock_t loop_s, loop_e;
 
-	// calculate the percentile (Matlab style)
+	// calculate the performance percentile (Matlab style)
 	qsort(q->tasks, q->N, sizeof(struct Task *), tasksort);
 	p = 0.95*q->N + 0.5;
 	pi = maximum(minimum(floor(p), q->N-1), 1);
@@ -111,7 +304,9 @@ void consistency_repeats(struct Queue *q, long repeats, TrainType traintype)
 	boundary = (1 - pr)*q->tasks[((long) pi)-1]->performance;
 	boundary += pr*q->tasks[((long) pi)]->performance;
 	note("boundary determined at: %f\n", boundary);
-
+	
+	// find the number of tasks that perform at least as good as the 95th
+	// percentile
 	N = 0;
 	for (i=0; i<q->N; i++)
 		if (q->tasks[i]->performance >= boundary)
@@ -121,12 +316,14 @@ void consistency_repeats(struct Queue *q, long repeats, TrainType traintype)
 	mean = Calloc(double, N);
 	perf = Calloc(double, N*repeats);
 
+	// create a new task queue with the tasks which perform well
 	nq->tasks = Malloc(struct Task *, N);
 	for (i=q->N-1; i>q->N-N-1; i--) 
 		nq->tasks[q->N-i-1] = q->tasks[i];
 	nq->N = N;
 	nq->i = 0;
 
+	// for each task run the consistency repeats
 	for (i=0; i<N; i++) {
 		task = get_next_task(nq);
 		make_model_from_task(task, model);
@@ -140,7 +337,8 @@ void consistency_repeats(struct Queue *q, long repeats, TrainType traintype)
 		for (r=0; r<repeats; r++) {
 			if (traintype == CV) {
 				loop_s = clock();
-				p = cross_validation(model, NULL, task->train_data, task->folds);
+				p = cross_validation(model, NULL, 
+						task->train_data, task->folds);
 				loop_e = clock();
 				time += elapsed_time(loop_s, loop_e);
 				matrix_set(perf, repeats, i, r, p);
@@ -152,15 +350,24 @@ void consistency_repeats(struct Queue *q, long repeats, TrainType traintype)
 			note("%3.3f\t", p);
 		}
 		for (r=0; r<repeats; r++) {
-			std[i] += pow(matrix_get(perf, repeats, i, r) - mean[i], 2);
+			std[i] += pow(matrix_get(
+						perf, 
+						repeats, 
+						i, 
+						r) - mean[i],
+				       	2.0);
 		}
 		std[i] /= ((double) repeats) - 1.0;
 		std[i] = sqrt(std[i]);
-		note("(m = %3.3f, s = %3.3f, t = %3.3f)\n", mean[i], std[i], time);
+		note("(m = %3.3f, s = %3.3f, t = %3.3f)\n", 
+				mean[i], std[i], time);
 	}
 
+	// find the best overall configurations: those with high average
+	// performance and low deviation in the performance
 	note("\nBest overall configuration(s):\n");
-	note("ID\tweights\tepsilon\t\tp\t\tkappa\t\tlambda\t\tmean_perf\tstd_perf\n");
+	note("ID\tweights\tepsilon\t\tp\t\tkappa\t\tlambda\t\t"
+			"mean_perf\tstd_perf\n");
 	p = 0.0;
 	bool breakout = false;
 	while (breakout == false) {
@@ -168,13 +375,17 @@ void consistency_repeats(struct Queue *q, long repeats, TrainType traintype)
 		pr = prctile(std, N, p/100.0);
 		for (i=0; i<N; i++) 
 			if ((pi - mean[i] < 0.0001) && (std[i] - pr < 0.0001)) {
-				note("(%li)\tw = %li\te = %f\tp = %f\tk = %f\tl = %f\t"
+				note("(%li)\tw = %li\te = %f\tp = %f\t"
+						"k = %f\tl = %f\t"
 						"mean: %3.3f\tstd: %3.3f\n",
 						nq->tasks[i]->ID, 
 						nq->tasks[i]->weight_idx,
-						nq->tasks[i]->epsilon, nq->tasks[i]->p,
-						nq->tasks[i]->kappa, nq->tasks[i]->lambda,
-						mean[i], std[i]);
+						nq->tasks[i]->epsilon, 
+						nq->tasks[i]->p,
+						nq->tasks[i]->kappa, 
+						nq->tasks[i]->lambda,
+						mean[i], 
+						std[i]);
 				breakout = true;
 			}
 		p += 1.0;
@@ -187,6 +398,30 @@ void consistency_repeats(struct Queue *q, long repeats, TrainType traintype)
 	free(mean);
 }
 
+/**
+ * @brief Run cross validation with a seed model
+ *
+ * @details
+ * This is an implementation of cross validation which uses the optimal
+ * parameters MajModel::V of a previous fold as initial conditions for
+ * MajModel::V of the next fold. An initial seed for V can be given through the
+ * seed_model parameter. If seed_model is NULL, random starting values are
+ * used.
+ *
+ * @todo
+ * The seed model shouldn't have to be allocated completely, since only V is
+ * used.
+ * @todo
+ * There must be some inefficiencies here because the fold model is allocated
+ * at every fold. This would be detrimental with large datasets.
+ *
+ * @param[in] 	model 		MajModel with the configuration to train
+ * @param[in] 	seed_model 	MajModel with a seed for MajModel::V
+ * @param[in] 	data 		MajData with the dataset
+ * @param[in] 	folds 		number of cross validation folds
+ * @returns 			performance (hitrate) of the configuration on 
+ * 				cross validation
+ */
 double cross_validation(struct MajModel *model, struct MajModel *seed_model,
 		struct MajData *data, long folds) 
 {
@@ -202,7 +437,7 @@ double cross_validation(struct MajModel *model, struct MajModel *seed_model,
 	double *performance = Calloc(double, folds);
 
 	if (seed_model == NULL) {
-		seed_model = Malloc(struct MajModel, 1);
+		seed_model = msvmmaj_init_model();
 		seed_model->n = 0; // we never use anything other than V
 		seed_model->m = model->m;
 		seed_model->K = model->K;
@@ -211,34 +446,40 @@ double cross_validation(struct MajModel *model, struct MajModel *seed_model,
 		fs = true;
 	}
 
-	train_data = Malloc(struct MajData, 1);
-	test_data = Malloc(struct MajData, 1);
-
+	train_data = msvmmaj_init_data();
+	test_data = msvmmaj_init_data();
+ 	// create splits
 	msvmmaj_make_cv_split(model->n, folds, cv_idx);
+
 	for (f=0; f<folds; f++) {
 		msvmmaj_get_tt_split(data, train_data, test_data, cv_idx, f);
-
-		fold_model = Malloc(struct MajModel, 1);
+		// initialize a model for this fold and copy the model
+		// parameters
+		fold_model = msvmmaj_init_model();
 		copy_model(model, fold_model);
 	
 		fold_model->n = train_data->n;
 		fold_model->m = train_data->m;
 		fold_model->K = train_data->K;
-		
+	
+		// allocate, initialize and seed the fold model	
 		msvmmaj_allocate_model(fold_model);
 		msvmmaj_initialize_weights(train_data, fold_model);
 		msvmmaj_seed_model_V(seed_model, fold_model);
-	
+
+		// train the model (without output)	
 		fid = MSVMMAJ_OUTPUT_FILE;
 		MSVMMAJ_OUTPUT_FILE = NULL;
 		msvmmaj_optimize(fold_model, train_data);
 		MSVMMAJ_OUTPUT_FILE = fid;
 
+		// calculate predictive performance on test set 
 		predy = Calloc(long, test_data->n);
 		msvmmaj_predict_labels(test_data, fold_model, predy);
 		performance[f] = msvmmaj_prediction_perf(test_data, predy);
 		total_perf += performance[f]/((double) folds);
 
+		// seed the seed model with the fold model
 		msvmmaj_seed_model_V(fold_model, seed_model);
 		
 		free(predy);
@@ -250,6 +491,7 @@ double cross_validation(struct MajModel *model, struct MajModel *seed_model,
 		msvmmaj_free_model(fold_model);
 	}
 
+	// if a seed model was allocated before, free it.
 	if (fs)
 		msvmmaj_free_model(seed_model);
 	free(train_data);
@@ -261,12 +503,28 @@ double cross_validation(struct MajModel *model, struct MajModel *seed_model,
 
 }
 
+/**
+ * @brief Run the grid search for a cross validation dataset
+ *
+ * @details
+ * Given a Queue of Task struct to be trained, a grid search is launched to
+ * find the optimal parameter configuration. As is also done within
+ * cross_validation(), the optimal weights of one parameter set are used as
+ * initial estimates for MajModel::V in the next parameter set. Note that to
+ * optimally exploit this feature of the optimization algorithm, the order in
+ * which tasks are considered is important. This is considered in 
+ * make_queue().
+ * 
+ * The performance found by cross validation is stored in the Task struct. 
+ *
+ * @param[in,out] 	q 	Queue with Task instances to run
+ */
 void start_training_cv(struct Queue *q)
 {
 	double perf, current_max = 0;
 	struct Task *task = get_next_task(q);
-	struct MajModel *seed_model = Malloc(struct MajModel, 1);
-	struct MajModel *model = Malloc(struct MajModel, 1);
+	struct MajModel *seed_model = msvmmaj_init_model();
+	struct MajModel *model = msvmmaj_init_model();
 	clock_t main_s, main_e, loop_s, loop_e;
 
 	model->n = task->train_data->n;
@@ -282,13 +540,16 @@ void start_training_cv(struct Queue *q)
 	
 	main_s = clock();
 	while (task) {
-		note("(%03li/%03li)\tw = %li\te = %f\tp = %f\tk = %f\t l = %f\t",
-				task->ID+1, q->N, task->weight_idx, task->epsilon,
+		note("(%03li/%03li)\tw = %li\te = %f\tp = %f\tk = %f\t "
+				"l = %f\t",
+				task->ID+1, q->N, task->weight_idx,
+			       	task->epsilon,
 				task->p, task->kappa, task->lambda);
 		make_model_from_task(task, model);
 		
 		loop_s = clock();
-		perf = cross_validation(model, seed_model, task->train_data, task->folds);
+		perf = cross_validation(model, seed_model, task->train_data,
+			       	task->folds);
 		loop_e = clock();
 		current_max = maximum(current_max, perf);
 
@@ -308,6 +569,23 @@ void start_training_cv(struct Queue *q)
 	msvmmaj_free_model(seed_model);
 }
 
+/**
+ * @brief Run the grid search for a train/test dataset
+ *
+ * @details
+ * This function is similar to start_training_cv(), except that the
+ * pre-determined training set is used only once, and the pre-determined test
+ * set is used for validation.
+ *
+ * @todo
+ * It would probably be better to train the model on the training set using
+ * cross validation and only use the test set when comparing with other
+ * methods. The way it is now, you're finding out which parameters predict
+ * _this_ test set best, which is not what you want.
+ *
+ * @param[in] 	q 	Queue with Task structs to run
+ *
+ */
 void start_training_tt(struct Queue *q)
 {
 	FILE *fid;
@@ -317,7 +595,7 @@ void start_training_tt(struct Queue *q)
 	double total_perf, current_max = 0;
 
 	struct Task *task = get_next_task(q);
-	struct MajModel *seed_model = Malloc(struct MajModel, 1);
+	struct MajModel *seed_model = msvmmaj_init_model();
 
 	clock_t main_s, main_e;
 	clock_t loop_s, loop_e;
@@ -334,7 +612,7 @@ void start_training_tt(struct Queue *q)
 				c+1, q->N, task->weight_idx, task->epsilon,
 				task->p, task->kappa, task->lambda);
 		loop_s = clock();
-		struct MajModel *model = Malloc(struct MajModel, 1);
+		struct MajModel *model = msvmmaj_init_model();
 		make_model_from_task(task, model);
 
 		model->n = task->train_data->n;
@@ -374,15 +652,37 @@ void start_training_tt(struct Queue *q)
 	msvmmaj_free_model(seed_model);
 }
 
+/**
+ * @brief Free the Queue struct
+ *
+ * @details
+ * Freeing the allocated memory of the Queue means freeing every Task struct
+ * and then freeing the Queue.
+ *
+ * @param[in] 	q 	Queue to be freed
+ *
+ */
 void free_queue(struct Queue *q)
 {
 	long i;
-	for (i=0; i<q->N; i++)
+	for (i=0; i<q->N; i++) {
+		free(q->tasks[i]->kernel_param);
 		free(q->tasks[i]);
+	}
 	free(q->tasks);
 	free(q);
 }
 
+/**
+ * @brief Copy parameters from Task to MajModel
+ *
+ * @details
+ * A Task struct only contains the parameters of the MajModel to be estimated.
+ * This function is used to copy these parameters.
+ *
+ * @param[in] 		task 	Task instance with parameters
+ * @param[in,out] 	model 	MajModel to which the parameters are copied
+ */
 void make_model_from_task(struct Task *task, struct MajModel *model)
 {
 	model->weight_idx = task->weight_idx;
@@ -392,6 +692,16 @@ void make_model_from_task(struct Task *task, struct MajModel *model)
 	model->lambda = task->lambda;
 }
 
+/**
+ * @brief Copy model parameters between two MajModel structs
+ *
+ * @details
+ * The parameters copied are MajModel::weight_idx, MajModel::epsilon,
+ * MajModel::p, MajModel::kappa, and MajModel::lambda.
+ *
+ * @param[in] 		from 	MajModel to copy parameters from
+ * @param[in,out] 	to 	MajModel to copy parameters to
+ */
 void copy_model(struct MajModel *from, struct MajModel *to)
 {
 	to->weight_idx = from->weight_idx;
diff --git a/src/mylapack.c b/src/mylapack.c
deleted file mode 100644
index 4a9cf81..0000000
--- a/src/mylapack.c
+++ /dev/null
@@ -1,49 +0,0 @@
-/**
- * @file mylapack.c
- * @author Gertjan van den Burg (burg@ese.eur.nl)
- * @date August 9, 2013
- * @brief Utility functions for interacting with LAPACK
- *
- * @details
- * Functions in this file are auxiliary functions which make it easier
- * to use LAPACK functions from liblapack.
- */
-
-#include "mylapack.h"
-
-/**
- * @name dposv
- * @brief Solve a system of equations AX = B where A is symmetric positive definite.
- * @ingroup libMSVMMaj
- *
- * See the LAPACK documentation at: 
- * http://www.netlib.org/lapack/explore-html/dc/de9/group__double_p_osolve.html
- */
-int dposv(char UPLO, int N, int NRHS, double *A, int LDA, double *B,
-		int LDB)
-{
-	extern void dposv_(char *UPLO, int *Np, int *NRHSp, double *A,
-			int *LDAp, double *B, int *LDBp, int *INFOp);
-	int INFO;
-	dposv_(&UPLO, &N, &NRHS, A, &LDA, B, &LDB, &INFO);
-	return INFO;
-}
-
-/**
- * @name dsysv
- * @brief Solve a system of equations AX = B where A is symmetric.
- * @ingroup libMSVMMaj
- *
- * See the LAPACK documentation at:
- * http://www.netlib.org/lapack/explore-html/d6/d0e/group__double_s_ysolve.html
- */
-int dsysv(char UPLO, int N, int NRHS, double *A, int LDA, int *IPIV,
-		double *B, int LDB, double *WORK, int LWORK)
-{
-	extern void dsysv_(char *UPLO, int *Np, int *NRHSp, double *A,
-			int *LDAp, int *IPIV, double *B, int *LDBp,
-			double *WORK, int *LWORK, int *INFOp);
-	int INFO;
-	dsysv_(&UPLO, &N, &NRHS, A, &LDA, IPIV, B, &LDB, WORK, &LWORK, &INFO);
-	return INFO;
-}
diff --git a/src/predMSVMMaj.c b/src/predMSVMMaj.c
index 966c7c0..3e3a101 100644
--- a/src/predMSVMMaj.c
+++ b/src/predMSVMMaj.c
@@ -1,17 +1,42 @@
+/**
+ * @file predMSVMMaj.c
+ * @author Gertjan van den Burg
+ * @date January, 2014
+ * @brief Command line interface for predicting class labels
+ *
+ * @details
+ * This is a command line program for predicting the class labels or
+ * determining the predictive performance of a pre-determined model on a given
+ * test dataset. The predictive performance can be written to the screen or
+ * the predicted class labels can be written to a specified output file. This
+ * is done using msvmmaj_write_predictions().
+ *
+ * The specified model file must follow the specification given in
+ * msvmmaj_write_model().
+ *
+ * For usage information, see the program help function.
+ *
+ */
+
+#include "msvmmaj.h"
+#include "msvmmaj_init.h"
 #include "msvmmaj_pred.h"
-#include "MSVMMaj.h"
 #include "util.h"
 
 #define MINARGS 3
 
 extern FILE *MSVMMAJ_OUTPUT_FILE;
 
+// function declarations
 void print_null(const char *s) {}
 void exit_with_help();
-void parse_command_line(int argc, char **argv, struct MajModel *model,
+void parse_command_line(int argc, char **argv, 
 		char *input_filename, char *output_filename, 
 		char *model_filename);
 
+/**
+ * @brief Help function
+ */
 void exit_with_help()
 {
 	printf("This is MSVMMaj, version %1.1f\n\n", VERSION);
@@ -22,6 +47,24 @@ void exit_with_help()
 	exit(0);
 }
 
+/**
+ * @brief Main interface function for predMSVMMaj
+ *
+ * @details
+ * Main interface for the command line program. A given model file is read and
+ * a test dataset is initialized from the given data. The predictive
+ * performance (hitrate) of the model on the test set is printed to the output
+ * stream (default = stdout). If an output file is specified the predictions
+ * are written to the file.
+ *
+ * @todo
+ * Ensure that the program can read model files without class labels
+ * specified. In that case no prediction accuracy is printed to the screen.
+ *
+ * @param[in] 	argc 	number of command line arguments
+ * @param[in] 	argv 	array of command line arguments
+ *
+ */
 int main(int argc, char **argv)
 {
 	long *predy;
@@ -31,16 +74,15 @@ int main(int argc, char **argv)
 	char model_filename[MAX_LINE_LENGTH];
 	char output_filename[MAX_LINE_LENGTH];;
 
-	struct MajModel *model = Malloc(struct MajModel, 1);
-	struct MajData *data = Malloc(struct MajData, 1);
-
 	if (argc < MINARGS || msvmmaj_check_argv(argc, argv, "-help") 
 			|| msvmmaj_check_argv_eq(argc, argv, "-h") )
 		exit_with_help();
-	parse_command_line(argc, argv, model, input_filename, output_filename,
+	parse_command_line(argc, argv, input_filename, output_filename,
 			model_filename);
 
-	// TODO: make sure that read_data allows for files without labels
+	// read the data and model
+	struct MajModel *model = msvmmaj_init_model();
+	struct MajData *data = msvmmaj_init_data();
 	msvmmaj_read_data(data, input_filename);
 	msvmmaj_read_model(model, model_filename);
 
@@ -50,8 +92,14 @@ int main(int argc, char **argv)
 				"does not equal the number of attributes in "
 				"model (%li)\n", data->m, model->m);
 		exit(1);
+	} else if (data->K != model->K) {
+		fprintf(stderr, "Error: number of classes in data (%li) "
+				"does not equal the number of classes in "
+				"model (%li)\n", data->K, model->K);
+		exit(1);
 	}
 
+	// predict labels and performance if test data has labels
 	predy = Calloc(long, data->n);
 	msvmmaj_predict_labels(data, model, predy);
 	if (data->y != NULL) {
@@ -59,11 +107,13 @@ int main(int argc, char **argv)
 		note("Predictive performance: %3.2f%%\n", performance);
 	}
 
+	// if output file is specified, write predictions to it
 	if (msvmmaj_check_argv_eq(argc, argv, "-o")) {
 		msvmmaj_write_predictions(data, predy, output_filename);
 		note("Predictions written to: %s\n", output_filename);
 	}
 
+	// free the model, data, and predictions
 	msvmmaj_free_model(model);
 	msvmmaj_free_data(data);
 	free(predy);
@@ -71,8 +121,26 @@ int main(int argc, char **argv)
 	return 0;
 }
 
-void parse_command_line(int argc, char **argv, struct MajModel *model, 
-		char *input_filename, char *output_filename, char *model_filename)
+/**
+ * @brief Parse command line arguments
+ *
+ * @details
+ * Read the data filename and model filename from the command line arguments.
+ * If specified, also read the output filename. If the quiet flag is given,
+ * set the global output stream to NULL. On error, exit_with_help().
+ *
+ * @param[in] 	argc 			number of command line arguments
+ * @param[in] 	argv 			array of command line arguments
+ * @param[in] 	input_filename 		pre-allocated array for the input 
+ * 					filename
+ * @param[in] 	output_filename 	pre-allocated array for the output
+ * 					filename
+ * @param[in] 	model_filename 		pre-allocated array for the model
+ * 					filename
+ *
+ */
+void parse_command_line(int argc, char **argv, char *input_filename,
+		char *output_filename, char *model_filename)
 {
 	int i;
 
@@ -91,7 +159,8 @@ void parse_command_line(int argc, char **argv, struct MajModel *model,
 				i--;
 				break;
 			default:
-				fprintf(stderr, "Unknown option: -%c\n", argv[i-1][1]);
+				fprintf(stderr, "Unknown option: -%c\n", 
+						argv[i-1][1]);
 				exit_with_help();
 		}
 	}
diff --git a/src/strutil.c b/src/strutil.c
index ae96239..ca4181f 100644
--- a/src/strutil.c
+++ b/src/strutil.c
@@ -1,5 +1,24 @@
+/**
+ * @file strutil.c
+ * @author Gertjan van den Burg
+ * @date January, 2014
+ * @brief Utility functions for dealing with strings
+ *
+ * @details
+ * This file contains functions for reading files, reading strings from a
+ * format and checking start and ends of strings.
+ */
+
 #include "strutil.h"
 
+/**
+ * @brief Check if a string starts with a prefix
+ *
+ * @param[in] 	str 	string
+ * @param[in] 	pre 	prefix
+ * @returns 		boolean, true if string starts with prefix, false
+ * 			otherwise
+ */
 bool str_startswith(const char *str, const char *pre)
 {
 	size_t lenpre = strlen(pre),
@@ -7,19 +26,41 @@ bool str_startswith(const char *str, const char *pre)
 	return lenstr < lenpre ? false : strncmp(pre, str, lenpre) == 0;
 }
 
+/**
+ * @brief Check if a string ends with a suffix
+ *
+ * @param[in] 	str 	string
+ * @param[in] 	suf 	suffix
+ * @returns 		boolean, true if string ends with suffix, false
+ * 			otherwise
+ */
 bool str_endswith(const char *str, const char *suf)
 {
 	size_t lensuf = strlen(suf),
 	       lenstr = strlen(str);
-	return lenstr < lensuf ? false : strncmp(str + lenstr - lensuf, suf, lensuf) == 0;
+	return lenstr < lensuf ? false : strncmp(str + lenstr - lensuf, suf, 
+			lensuf) == 0;
 }
 
+/**
+ * @brief Move to next line in file
+ *
+ * @param[in] 	fid 		File opened for reading
+ * @param[in] 	filename 	name of the file pointed to by fid
+ */
 void next_line(FILE *fid, char *filename)
 {
 	char buffer[MAX_LINE_LENGTH];
 	get_line(fid, filename, buffer);
 }
 
+/**
+ * @brief Read line to buffer
+ *
+ * @param[in] 		fid 		File opened for reading
+ * @param[in] 		filename 	name of the file
+ * @param[in,out] 	buffer 		allocated buffer to read to
+ */
 void get_line(FILE *fid, char *filename, char *buffer)
 {
 	if (fgets(buffer, MAX_LINE_LENGTH, fid) == NULL) {
@@ -28,6 +69,14 @@ void get_line(FILE *fid, char *filename, char *buffer)
 	}
 }
 
+/**
+ * @brief Read a double from file following a format
+ *
+ * @param[in] 	fid 		File opened for reading
+ * @param[in] 	filename 	Name of the file
+ * @param[in] 	fmt 		Format containing a float format
+ * @returns 			value read (if any)
+ */
 double get_fmt_double(FILE *fid, char *filename, const char *fmt)
 {
 	char buffer[MAX_LINE_LENGTH];
@@ -38,6 +87,14 @@ double get_fmt_double(FILE *fid, char *filename, const char *fmt)
 	return value;
 }
 
+/**
+ * @brief Read a long integer from file following a format
+ *
+ * @param[in] 	fid 		File opened for reading
+ * @param[in] 	filename 	Name of the file
+ * @param[in] 	fmt 		Format containing a long integer format
+ * @returns 			value read (if any)
+ */
 long get_fmt_long(FILE *fid, char *filename, const char *fmt)
 {
 	char buffer[MAX_LINE_LENGTH];
@@ -48,6 +105,20 @@ long get_fmt_long(FILE *fid, char *filename, const char *fmt)
 	return value;
 }
 
+/**
+ * @brief Read all doubles in a given buffer
+ *
+ * @details
+ * This function is used to read a line of doubles from a buffer. All the
+ * doubles found are stored in a pre-allocated array. 
+ *
+ * @param[in] 	buffer 		a string buffer
+ * @param[in] 	offset 		an offset of the string to start looking for
+ * 				doubles
+ * @param[in] 	all_doubles 	pre-allocated array of doubles (should be large
+ * 				enough)
+ * @returns 			number of doubles read
+ */
 long all_doubles_str(char *buffer, long offset, double *all_doubles)
 {
 	double value;
@@ -69,6 +140,20 @@ long all_doubles_str(char *buffer, long offset, double *all_doubles)
 	return i;
 }
 
+/**
+ * @brief Read all longs in a given buffer
+ *
+ * @details
+ * This function is used to read a line of longs from a buffer. All the
+ * longs found are stored in a pre-allocated array. 
+ *
+ * @param[in] 	buffer 		a string buffer
+ * @param[in] 	offset 		an offset of the string to start looking for
+ * 				longs
+ * @param[in] 	all_longs 	pre-allocated array of longs (should be large
+ * 				enough)
+ * @returns 			number of longs read
+ */
 long all_longs_str(char *buffer, long offset, long *all_longs)
 {
 	long value;
diff --git a/src/timer.c b/src/timer.c
index 2187fb2..3a763a0 100644
--- a/src/timer.c
+++ b/src/timer.c
@@ -1,7 +1,25 @@
+/**
+ * @file timer.c
+ * @author Gertjan van den Burg
+ * @date January, 2014
+ * @brief Function for calculating time difference
+ *
+ * @details
+ * This file contains a simple function for calculating the time in seconds
+ * elapsed between two clock() calls.
+ */
+
 #include <time.h>
 
 #include "timer.h"
 
+/**
+ * @brief Calculate the time between two clocks
+ *
+ * @param[in] 	s_time 	starting time
+ * @param[in] 	e_time 	end time
+ * @returns 		time elapsed in seconds
+ */
 double elapsed_time(clock_t s_time, clock_t e_time)
 {
 	return ((double) (e_time - s_time))/((double) CLOCKS_PER_SEC);
diff --git a/src/trainMSVMMaj.c b/src/trainMSVMMaj.c
index b4b74df..e045a6c 100644
--- a/src/trainMSVMMaj.c
+++ b/src/trainMSVMMaj.c
@@ -1,54 +1,93 @@
+/**
+ * @file trainMSVMMaj.c
+ * @author Gertjan van den Burg
+ * @date August, 2013
+ * @brief Command line interface for training a single model with MSVMMaj
+ *
+ * @details
+ * This is a command line program for training a single model on a given
+ * dataset. To run a grid search over a number of parameter configurations, 
+ * see trainMSVMMajdataset.c. 
+ *
+ */
+
 #include <time.h>
 #include <math.h>
 
+#include "msvmmaj_kernel.h"
 #include "libMSVMMaj.h"
+#include "msvmmaj.h"
+#include "msvmmaj_init.h"
 #include "msvmmaj_train.h"
 #include "util.h"
-#include "MSVMMaj.h"
 
 #define MINARGS 2
 
 extern FILE *MSVMMAJ_OUTPUT_FILE;
 
+// function declarations
 void print_null(const char *s) {}
 void exit_with_help();
 void parse_command_line(int argc, char **argv, struct MajModel *model, 
 		char *input_filename, char *output_filename, char *model_filename);
 
+/**
+ * @brief Help function
+ */
 void exit_with_help()
 {
 	printf("This is MSVMMaj, version %1.1f\n\n", VERSION);
 	printf("Usage: trainMSVMMaj [options] training_data_file\n");
 	printf("Options:\n");
+	printf("-c coef : coefficient for the polynomial and sigmoid kernel\n");
+	printf("-d degree : degree for the polynomial kernel\n");
 	printf("-e epsilon : set the value of the stopping criterion\n");
+	printf("-g gamma : parameter for the rbf, polynomial or sigmoid "
+			"kernel\n");
 	printf("-h | -help : print this help.\n");
 	printf("-k kappa : set the value of kappa used in the Huber hinge\n");
 	printf("-l lambda : set the value of lambda (lambda > 0)\n");
 	printf("-m model_file : use previous model as seed for W and t\n");
 	printf("-o output_file : write output to file\n");
-	printf("-p p-value : set the value of p in the lp norm (1.0 <= p <= 2.0)\n");
+	printf("-p p-value : set the value of p in the lp norm "
+			"(1.0 <= p <= 2.0)\n");
 	printf("-q : quiet mode (no output)\n");
-	printf("-r rho : choose the weigth specification (1 = unit, 2 = group)\n");
+	printf("-r rho : choose the weigth specification (1 = unit, 2 = "
+			"group)\n");
+	printf("-t type: kerneltype (LINEAR=0, POLY=1, RBF=2, SIGMOID=3)\n");
+	printf("-u use_cholesky: use cholesky decomposition when using "
+			"kernels (0 = false, 1 = true). Default 0.\n");
 
 	exit(0);
 }
 
-/*
-	Main
-*/
+/**
+ * @brief Main interface function for trainMSVMMaj
+ *
+ * @details
+ * Main interface for the command line program. A given dataset file is read
+ * and a MSVMMaj model is trained on this data. By default the progress of the
+ * computations are written to stdout. See for full options of the program the
+ * help function.
+ *
+ * @param[in] 	argc 	number of command line arguments
+ * @param[in] 	argv 	array of command line arguments
+ *
+ */
 int main(int argc, char **argv)
 {
 	char input_filename[MAX_LINE_LENGTH];
 	char model_filename[MAX_LINE_LENGTH];
 	char output_filename[MAX_LINE_LENGTH];
 
-	struct MajModel *model = Malloc(struct MajModel, 1);
-	struct MajData *data = Malloc(struct MajData, 1);
+	struct MajModel *model = msvmmaj_init_model();
+	struct MajData *data = msvmmaj_init_data();
 
 	if (argc < MINARGS || msvmmaj_check_argv(argc, argv, "-help") 
 			|| msvmmaj_check_argv_eq(argc, argv, "-h") ) 
 		exit_with_help();
-	parse_command_line(argc, argv, model, input_filename, output_filename, model_filename);
+	parse_command_line(argc, argv, model, input_filename, 
+			output_filename, model_filename);
 
 	// read data file
 	msvmmaj_read_data(data, input_filename);
@@ -59,22 +98,25 @@ int main(int argc, char **argv)
 	model->K = data->K;
 	model->data_file = input_filename;
 
+	// initialize kernel (if necessary)
+	msvmmaj_make_kernel(model, data);
+
 	// allocate model and initialize weights
 	msvmmaj_allocate_model(model);
 	msvmmaj_initialize_weights(data, model);
 
+	// seed the random number generator (only place in programs is in
+	// command line interfaces)
 	srand(time(NULL));
 
 	if (msvmmaj_check_argv_eq(argc, argv, "-m")) {
-		struct MajModel *seed_model = Malloc(struct MajModel, 1);
+		struct MajModel *seed_model = msvmmaj_init_model();
 		msvmmaj_read_model(seed_model, model_filename);
 		msvmmaj_seed_model_V(seed_model, model);
 		msvmmaj_free_model(seed_model);
 	} else {
 		msvmmaj_seed_model_V(NULL, model);
 	}
-	// initialize kernel (if necessary)
-	// msvmmaj_make_kernel(model, data);
 
 	// start training
 	msvmmaj_optimize(model, data);
@@ -92,18 +134,34 @@ int main(int argc, char **argv)
 	return 0;
 }
 
+/**
+ * @brief Parse command line arguments
+ *
+ * @details
+ * Process the command line arguments for the model parameters, and record
+ * them in the specified MajModel. An input filename for the dataset is read
+ * and if specified an output filename and a model filename for the seed
+ * model.
+ *
+ * @param[in] 	argc 			number of command line arguments
+ * @param[in] 	argv 			array of command line arguments
+ * @param[in] 	model 			initialized model
+ * @param[in] 	input_filename 		pre-allocated buffer for the input
+ * 					filename
+ * @param[in] 	output_filename 	pre-allocated buffer for the output
+ * 					filename
+ * @param[in] 	model_filename 		pre-allocated buffer for the model
+ * 					filename
+ *
+ */
 void parse_command_line(int argc, char **argv, struct MajModel *model, 
 		char *input_filename, char *output_filename, char *model_filename)
 {
-	int i;
-
-	// default values
-	model->p = 1.0;
-	model->lambda = pow(2, -8.0);
-	model->epsilon = 1e-6;
-	model->kappa = 0.0;
-	model->weight_idx = 1;
-	
+	int i, tmp;
+	double gamma = 1.0, 
+	       degree = 2.0, 
+	       coef = 0.0;
+
 	MSVMMAJ_OUTPUT_FILE = stdout;
 
 	// parse options
@@ -113,9 +171,18 @@ void parse_command_line(int argc, char **argv, struct MajModel *model,
 			exit_with_help();
 		}
 		switch (argv[i-1][1]) {
+			case 'c':
+				coef = atof(argv[i]);
+				break;
+			case 'd':
+				degree = atof(argv[i]);
+				break;
 			case 'e':
 				model->epsilon = atof(argv[i]);
 				break;
+			case 'g':
+				gamma = atof(argv[i]);
+				break;
 			case 'k':
 				model->kappa = atof(argv[i]);
 				break;
@@ -134,20 +201,50 @@ void parse_command_line(int argc, char **argv, struct MajModel *model,
 			case 'r':
 				model->weight_idx = atoi(argv[i]);
 				break;
+			case 't':
+				model->kerneltype = atoi(argv[i]);
+				break;
+			case 'u':
+				tmp = atoi(argv[i]);
+				if (!(tmp == 1 || tmp == 0)) 
+					fprintf(stderr, "Unknown value %i for"
+							" use_cholesky", tmp);
+				model->use_cholesky = (tmp == 1) ? true : false;
+				break;
 			case 'q':
 				MSVMMAJ_OUTPUT_FILE = NULL;
 				i--;
 				break;
 			default:
-				fprintf(stderr, "Unknown option: -%c\n", argv[i-1][1]);
+				fprintf(stderr, "Unknown option: -%c\n", 
+						argv[i-1][1]);
 				exit_with_help();
 		}
 	}
-
+	
 	// read input filename
 	if (i >= argc)
 		exit_with_help();
 
 	strcpy(input_filename, argv[i]);
-}
 
+	// set kernel parameters
+	switch (model->kerneltype) {
+		case K_LINEAR:
+			break;
+		case K_POLY:
+			model->kernelparam = Calloc(double, 3);
+			model->kernelparam[0] = gamma;
+			model->kernelparam[1] = coef;
+			model->kernelparam[2] = degree;
+			break;
+		case K_RBF:
+			model->kernelparam = Calloc(double, 1);
+			model->kernelparam[0] = gamma;
+			break;
+		case K_SIGMOID:
+			model->kernelparam = Calloc(double, 1);
+			model->kernelparam[0] = gamma;
+			model->kernelparam[1] = coef;
+	}
+}
diff --git a/src/trainMSVMMajdataset.c b/src/trainMSVMMajdataset.c
index 7c3385c..097df85 100644
--- a/src/trainMSVMMajdataset.c
+++ b/src/trainMSVMMajdataset.c
@@ -1,7 +1,28 @@
+/**
+ * @file trainMSVMMajdataset.c
+ * @author Gertjan van den Burg
+ * @date January, 2014
+ * @brief Command line interface for the grid search program
+ *
+ * @details
+ * This is a command line interface to the parameter grid search functionality
+ * of the algorithm. The grid search is specified in a separate file, thereby
+ * reducing the number of command line arguments. See
+ * read_training_from_file() for documentation on the training file. 
+ *
+ * The program runs a grid search as specified in the training file. If
+ * desired the grid search can incorporate consistency checks to find the
+ * configuration among the best configurations which scores consistently high.
+ * All output is written to stdout, unless the quiet mode is specified.
+ *
+ * For further usage information, see the program help function.
+ *
+ */
+
 #include <time.h>
 
 #include "crossval.h"
-#include "MSVMMaj.h"
+#include "msvmmaj.h"
 #include "msvmmaj_pred.h"
 #include "msvmmaj_train.h"
 #include "msvmmaj_train_dataset.h"
@@ -12,11 +33,15 @@
 
 extern FILE *MSVMMAJ_OUTPUT_FILE;
 
+// function declarations
 void print_null(const char *s) {}
 void exit_with_help();
 void parse_command_line(int argc, char **argv, char *input_filename);
 void read_training_from_file(char *input_filename, struct Training *training);
 
+/**
+ * @brief Help function
+ */
 void exit_with_help()
 {
 	printf("This is MSVMMaj, version %1.1f\n\n", VERSION);
@@ -28,6 +53,22 @@ void exit_with_help()
 	exit(0);
 }
 
+/**
+ * @brief Main interface function for trainMSVMMajdataset
+ *
+ * @details
+ * Main interface for the command line program. A given training file which
+ * specifies a grid search over a single dataset is read. From this, a Queue
+ * is created containing all Task instances that need to be performed in the
+ * search. Depending on the type of dataset, either cross validation or
+ * train/test split training is performed for all tasks. If specified,
+ * consistency repeats are done at the end of the grid search. Note that
+ * currently no output is produced other than what is written to stdout.
+ *
+ * @param[in] 	argc 	number of command line arguments
+ * @param[in] 	argv 	array of command line arguments
+ *
+ */
 int main(int argc, char **argv)
 {
 	char input_filename[MAX_LINE_LENGTH];
@@ -78,6 +119,21 @@ int main(int argc, char **argv)
 	return 0;
 }
 
+/**
+ * @brief Parse command line arguments
+ *
+ * @details
+ * Few arguments can be supplied to the command line. Only quiet mode can be
+ * specified, or help can be requested. The filename of the training file is
+ * read from the arguments. Parsing of the training file is done separately in
+ * read_training_from_file().
+ *
+ * @param[in] 	argc 		number of command line arguments
+ * @param[in] 	argv 		array of command line arguments
+ * @param[in] 	input_filename 	pre-allocated buffer for the training
+ * 				filename.
+ *
+ */
 void parse_command_line(int argc, char **argv, char *input_filename)
 {
 	int i;
@@ -94,7 +150,8 @@ void parse_command_line(int argc, char **argv, char *input_filename)
 				i--;
 				break;
 			default:
-				fprintf(stderr, "Unknown option: -%c\n", argv[i-1][1]);
+				fprintf(stderr, "Unknown option: -%c\n", 
+						argv[i-1][1]);
 				exit_with_help();
 		}
 	}
@@ -105,6 +162,21 @@ void parse_command_line(int argc, char **argv, char *input_filename)
 	strcpy(input_filename, argv[i]);
 }
 
+/**
+ * @brief Read the Training struct from file
+ *
+ * @details
+ * Read the Training struct from a file. The training file follows a specific
+ * format specified in @ref spec_training_file. 
+ *
+ * Commonly used string functions in this function are all_doubles_str() and
+ * all_longs_str(). 
+ *
+ * @param[in] 	input_filename 	filename of the training file
+ * @param[in] 	training 	Training structure to place the parsed
+ * 				parameter grid.
+ *
+ */
 void read_training_from_file(char *input_filename, struct Training *training)
 {
 	long i, nr = 0;
@@ -117,7 +189,8 @@ void read_training_from_file(char *input_filename, struct Training *training)
 
 	fid = fopen(input_filename, "r");
 	if (fid == NULL) {
-		fprintf(stderr, "Error opening training file %s\n", input_filename);
+		fprintf(stderr, "Error opening training file %s\n", 
+				input_filename);
 		exit(1);
 	}
 	training->traintype = CV;
@@ -126,11 +199,13 @@ void read_training_from_file(char *input_filename, struct Training *training)
 		Memset(lparams, long, MAX_LINE_LENGTH);
 		if (str_startswith(buffer, "train:")) {
 			sscanf(buffer, "train: %s\n", train_filename);
-			training->train_data_file = Calloc(char, MAX_LINE_LENGTH);
+			training->train_data_file = Calloc(char, 
+					MAX_LINE_LENGTH);
 			strcpy(training->train_data_file, train_filename);
 		} else if (str_startswith(buffer, "test:")) {
 			sscanf(buffer, "test: %s\n", test_filename);
-			training->test_data_file = Calloc(char, MAX_LINE_LENGTH);
+			training->test_data_file = Calloc(char,
+				       	MAX_LINE_LENGTH);
 			strcpy(training->test_data_file, test_filename);
 			training->traintype = TT;
 		} else if (str_startswith(buffer, "p:")) {
@@ -167,16 +242,76 @@ void read_training_from_file(char *input_filename, struct Training *training)
 			nr = all_longs_str(buffer, 6, lparams);
 			training->folds = lparams[0];
 			if (nr > 1)
-				fprintf(stderr, "Field \"folds\" only takes one value. "
-						"Additional fields are ignored.\n");
+				fprintf(stderr, "Field \"folds\" only takes "
+						"one value. Additional "
+						"fields are ignored.\n");
 		} else if (str_startswith(buffer, "repeats:")) {
 			nr = all_longs_str(buffer, 8, lparams);
 			training->repeats = lparams[0];
 			if (nr > 1)
-				fprintf(stderr, "Field \"repeats\" only takes one value. "
-						"Additional fields are ignored.\n");
+				fprintf(stderr, "Field \"repeats\" only "
+						"takes one value. Additional "
+						"fields are ignored.\n");
+		} else if (str_startswith(buffer, "kernel:")) {
+			nr = all_longs_str(buffer, 7, lparams);
+			if (nr > 1)
+				fprintf(stderr, "Field \"kernel\" only takes "
+						"one value. Additional "
+						"fields are ignored.\n");
+			switch (lparams[0]) {
+				case 0:
+					training->kerneltype = K_LINEAR;
+					break;
+				case 1:
+					training->kerneltype = K_POLY;
+					break;
+				case 2:
+					training->kerneltype = K_RBF;
+					break;
+				case 3:
+					training->kerneltype = K_SIGMOID;
+					break;
+			}
+		} else if (str_startswith(buffer, "gamma:")) {
+			nr = all_doubles_str(buffer, 6, params);
+			if (training->kerneltype == K_LINEAR) {
+				fprintf(stderr, "Field \"gamma\" ignored, "
+						"linear kernel is used.\n");
+				training->Ng = 0;
+				break;
+			}
+			training->gammas = Calloc(double, nr);
+			for (i=0; i<nr; i++)
+				training->gammas[i] = params[i];
+			training->Ng = nr;
+		} else if (str_startswith(buffer, "coef:")) {
+			nr = all_doubles_str(buffer, 5, params);
+			if (training->kerneltype == K_LINEAR ||
+				training->kerneltype == K_RBF) {
+				fprintf(stderr, "Field \"coef\" ignored with"
+						"specified kernel.\n");
+				training->Nc = 0;
+				break;
+			}
+			training->coefs = Calloc(double, nr);
+			for (i=0; i<nr; i++)
+				training->coefs[i] = params[i];
+			training->Nc = nr;
+		} else if (str_startswith(buffer, "degree:")) {
+			nr = all_doubles_str(buffer, 7, params);
+			if (training->kerneltype != K_POLY) {
+				fprintf(stderr, "Field \"degree\" ignored "
+						"with specified kernel.\n");
+				training->Nd = 0;
+				break;
+			}
+			training->degrees = Calloc(double, nr);
+			for (i=0; i<nr; i++)
+				training->degrees[i] = params[i];
+			training->Nd = nr;
 		} else {
-			fprintf(stderr, "Cannot find any parameters on line: %s\n", buffer);
+			fprintf(stderr, "Cannot find any parameters on line: "
+					"%s\n", buffer);
 		}
 	}
 
diff --git a/src/util.c b/src/util.c
index 101cb00..8e4b806 100644
--- a/src/util.c
+++ b/src/util.c
@@ -1,19 +1,55 @@
+/**
+ * @file util.c
+ * @author Gertjan van den Burg
+ * @date January, 2014
+ * @brief Utility functions
+ *
+ * @details
+ * This file contains several utility functions for coordinating input and
+ * output of data and model files. It also contains string functions.
+ *
+ * @todo
+ * Pull this apart.
+ *
+ */
 #include <math.h>
 #include <stdarg.h>
 #include <time.h>
 
-#include "matrix.h"
-#include "MSVMMaj.h"
+#include "msvmmaj.h"
+#include "msvmmaj_matrix.h"
 #include "strutil.h"
 
 #include "util.h"
 
-FILE *MSVMMAJ_OUTPUT_FILE;
-
-/*
-	Read the data from the data_file. The data matrix X is augmented
-	with a column of ones, to get the matrix Z.
-*/
+FILE *MSVMMAJ_OUTPUT_FILE; 	///< The #MSVMMAJ_OUTPUT_FILE specifies the 
+				///< output stream to which all output is 
+				///< written. This is done through the
+				///< internal (!)
+				///< function msvmmaj_print_string(). The
+				///< advantage of using a global output 
+				///< stream variable is that the output can
+				///< temporarily be suppressed by importing
+				///< this variable through @c extern and 
+				///< (temporarily) setting it to NULL.
+
+/**
+ * @brief Read data from file
+ * 
+ * @details
+ * Read the data from the data_file. The data matrix X is augmented
+ * with a column of ones, to get the matrix Z. The data is expected
+ * to follow a specific format, which is specified in the @ref spec_data_file.
+ * The class labels are corrected internally to correspond to the interval 
+ * [1 .. K], where K is the total number of classes.
+ *
+ * @todo
+ * Make sure that this function allows datasets without class labels for
+ * testing.
+ *
+ * @param[in,out] 	dataset 	initialized MajData struct
+ * @param[in] 		data_file 	filename of the data file.
+ */
 void msvmmaj_read_data(struct MajData *dataset, char *data_file)
 {
 	FILE *fid;
@@ -22,7 +58,7 @@ void msvmmaj_read_data(struct MajData *dataset, char *data_file)
 	long nr = 0; // used to check consistency of data
 	double value;
 	long K = 0;
-	long min_y = 1000;
+	long min_y = 1000000;
 
 	char buf[MAX_LINE_LENGTH];
 
@@ -79,13 +115,15 @@ void msvmmaj_read_data(struct MajData *dataset, char *data_file)
 			dataset->y[i]++;
 		K++;
 	} else if (min_y < 0 ) {
-		fprintf(stderr, "ERROR: wrong class labels in %s, minimum value is: %ld\n",
+		fprintf(stderr, "ERROR: wrong class labels in %s, minimum "
+				"value is: %ld\n",
 				data_file, min_y);
 		exit(0);
 	}
 
 	if (nr < n * m) {
-		fprintf(stderr, "ERROR: not enough data found in %s\n", data_file);
+		fprintf(stderr, "ERROR: not enough data found in %s\n", 
+				data_file);
 		exit(0);
 	}
 	
@@ -98,6 +136,19 @@ void msvmmaj_read_data(struct MajData *dataset, char *data_file)
 	dataset->K = K;
 }
 
+/**
+ * @brief Read model from file
+ *
+ * @details
+ * Read a MajModel from a model file. The MajModel struct must have been
+ * initalized elswhere. The model file is expected to follow the @ref
+ * spec_model_file. The easiest way to generate a model file is through
+ * msvmmaj_write_model(), which can for instance be used in trainMSVMMaj.c.
+ *
+ * @param[in,out] 	model 		initialized MajModel
+ * @param[in] 		model_filename 	filename of the model file
+ *
+ */
 void msvmmaj_read_model(struct MajModel *model, char *model_filename)
 {
 	long i, j, nr = 0;
@@ -108,7 +159,8 @@ void msvmmaj_read_model(struct MajModel *model, char *model_filename)
 
 	fid = fopen(model_filename, "r");
 	if (fid == NULL) {
-		fprintf(stderr, "Error opening model file %s\n", model_filename);
+		fprintf(stderr, "Error opening model file %s\n", 
+				model_filename);
 		exit(1);
 	}
 	// skip the first four lines
@@ -120,7 +172,8 @@ void msvmmaj_read_model(struct MajModel *model, char *model_filename)
 	model->lambda = get_fmt_double(fid, model_filename, "lambda = %lf");
 	model->kappa = get_fmt_double(fid, model_filename, "kappa = %lf");
 	model->epsilon = get_fmt_double(fid, model_filename, "epsilon = %lf");
-	model->weight_idx = (int) get_fmt_long(fid, model_filename, "weight_idx = %li");
+	model->weight_idx = (int) get_fmt_long(fid, model_filename, 
+			"weight_idx = %li");
 
 	// skip to data section
 	for (i=0; i<2; i++)
@@ -128,7 +181,8 @@ void msvmmaj_read_model(struct MajModel *model, char *model_filename)
 
 	// read filename of data file
 	if (fgets(buffer, MAX_LINE_LENGTH, fid) == NULL) {
-		fprintf(stderr, "Error reading model file %s\n", model_filename);
+		fprintf(stderr, "Error reading model file %s\n", 
+				model_filename);
 		exit(1);
 	}
 	sscanf(buffer, "filename = %s\n", data_filename);
@@ -153,12 +207,25 @@ void msvmmaj_read_model(struct MajModel *model, char *model_filename)
 	}
 	if (nr != (model->m+1)*(model->K-1)) {
 		fprintf(stderr, "Error reading model file %s. "
-				"Not enough elements of V found.\n", model_filename);
+				"Not enough elements of V found.\n", 
+				model_filename);
 		exit(1);
 	}
-
 }
 
+/**
+ * @brief Write model to file
+ *
+ * @details
+ * Write a MajModel to a file. The current time is specified in the file in
+ * UTC + offset. The model file further corresponds to the @ref
+ * spec_model_file.
+ *
+ * @param[in] 	model 		MajModel which contains an estimate for 
+ * 				MajModel::V
+ * @param[in] 	output_filename the output file to write the model to
+ *
+ */
 void msvmmaj_write_model(struct MajModel *model, char *output_filename)
 {
 	FILE *fid;
@@ -171,7 +238,8 @@ void msvmmaj_write_model(struct MajModel *model, char *output_filename)
 	// open output file
 	fid = fopen(output_filename, "w");
 	if (fid == NULL) {
-		fprintf(stderr, "Error opening output file %s", output_filename);
+		fprintf(stderr, "Error opening output file %s", 
+				output_filename);
 		exit(1);
 	}
 
@@ -201,7 +269,8 @@ void msvmmaj_write_model(struct MajModel *model, char *output_filename)
 
 	// Write output to file
 	fprintf(fid, "Output file for MSVMMaj (version %1.1f)\n", VERSION);
-	fprintf(fid, "Generated on: %s (UTC %+03i:%02i)\n\n", timestr, hours, minutes);
+	fprintf(fid, "Generated on: %s (UTC %+03i:%02i)\n\n", 
+			timestr, hours, minutes);
 	fprintf(fid, "Model:\n");
 	fprintf(fid, "p = %15.16f\n", model->p);
 	fprintf(fid, "lambda = %15.16f\n", model->lambda);
@@ -218,35 +287,71 @@ void msvmmaj_write_model(struct MajModel *model, char *output_filename)
 	fprintf(fid, "Output:\n");
 	for (i=0; i<model->m+1; i++) {
 		for (j=0; j<model->K-1; j++) {
-			fprintf(fid, "%+15.16f ", matrix_get(model->V, model->K-1, i, j));
+			fprintf(fid, "%+15.16f ", 
+					matrix_get(model->V, 
+						model->K-1, i, j));
 		}
 		fprintf(fid, "\n");
 	}
 
 	fclose(fid);
-
 }
 
-void msvmmaj_write_predictions(struct MajData *data, long *predy, char *output_filename)
+/**
+ * @brief Write predictions to file
+ *
+ * @details
+ * Write the given predictions to an output file, such that the resulting file
+ * corresponds to the @ref spec_data_file. 
+ *
+ * @param[in] 	data 		MajData with the original instances
+ * @param[in] 	predy 		predictions of the class labels of the 
+ * 				instances in the given MajData. Note that the
+ * 				order of the instances is assumed to be the
+ * 				same.
+ * @param[in] 	output_filename the file to which the predictions are written
+ *
+ */
+void msvmmaj_write_predictions(struct MajData *data, long *predy, 
+		char *output_filename)
 {
 	long i, j;
 	FILE *fid;
 
 	fid = fopen(output_filename, "w");
 	if (fid == NULL) {
-		fprintf(stderr, "Error opening output file %s", output_filename);
+		fprintf(stderr, "Error opening output file %s", 
+				output_filename);
 		exit(1);
 	}
 
 	for (i=0; i<data->n; i++) {
 		for (j=0; j<data->m; j++) 
-			fprintf(fid, "%f ", matrix_get(data->Z, data->m+1, i, j+1));
+			fprintf(fid, "%f ", 
+					matrix_get(data->Z, 
+						data->m+1, i, j+1));
 		fprintf(fid, "%li\n", predy[i]);
 	}
 
 	fclose(fid);
 }
 
+/**
+ * @brief Check if any command line arguments contain string
+ *
+ * @details
+ * Check if any of a given array of command line arguments contains a given 
+ * string. If the string is found, the index of the string in argv is 
+ * returned. If the string is not found, 0 is returned.
+ *
+ * This function is copied from MSVMpack/libMSVM.c.
+ *
+ * @param[in] 	argc 	number of command line arguments
+ * @param[in] 	argv 	command line arguments
+ * @param[in] 	str 	string to find in the arguments
+ * @returns 		index of the string in the arguments if found, 0
+ * 			otherwise
+ */
 int msvmmaj_check_argv(int argc, char **argv, char *str)
 {
 	int i;
@@ -260,6 +365,22 @@ int msvmmaj_check_argv(int argc, char **argv, char *str)
 	return arg_str;
 }
 
+/**
+ * @brief Check if a command line argument equals a string
+ *
+ * @details
+ * Check if any of the command line arguments is exactly equal to a given
+ * string. If so, return the index of the corresponding command line argument.
+ * If not, return 0.
+ *
+ * This function is copied from MSVMpack/libMSVM.c
+ *
+ * @param[in] 	argc 	number of command line arguments
+ * @param[in] 	argv 	command line arguments
+ * @param[in] 	str 	string to find in the arguments
+ * @returns 		index of the command line argument that corresponds to
+ * 			the string, 0 if none matches.
+ */
 int msvmmaj_check_argv_eq(int argc, char **argv, char *str) 
 {
 	int i;
@@ -274,6 +395,19 @@ int msvmmaj_check_argv_eq(int argc, char **argv, char *str)
 }
 
 
+/**
+ * @brief Print a given string to the specified output stream
+ *
+ * @details
+ * This function is used to print a given string to the output stream
+ * specified by #MSVMMAJ_OUTPUT_FILE. The stream is flushed after the string 
+ * is written to the stream. If #MSVMMAJ_OUTPUT_FILE is NULL, nothing is
+ * written. Note that this function is only used by note(), it should never be
+ * used directly.
+ *
+ * @param[in] 	s 	string to write to the stream
+ *
+ */
 static void msvmmaj_print_string(const char *s)
 {
 	if (MSVMMAJ_OUTPUT_FILE != NULL) {
@@ -282,6 +416,19 @@ static void msvmmaj_print_string(const char *s)
 	}
 }
 
+/**
+ * @brief Parse a formatted string and write to the output stream
+ *
+ * @details
+ * This function is a replacement of fprintf(), such that the output stream
+ * does not have to be specified at each function call. The functionality is
+ * exactly the same however. Writing the formatted string to the output stream
+ * is handled by msvmmaj_print_string().
+ *
+ * @param[in] 	fmt 	String format
+ * @param[in] 	... 	variable argument list for the string format
+ *
+ */
 void note(const char *fmt,...)
 {
 	char buf[BUFSIZ];
@@ -292,6 +439,16 @@ void note(const char *fmt,...)
 	(*msvmmaj_print_string)(buf);
 }
 
+/**
+ * @brief Allocate memory for a MajModel
+ *
+ * @details
+ * This function can be used to allocate the memory needed for a MajModel. All
+ * arrays in the model are specified and initialized to 0.
+ *
+ * @param[in] 	model 	MajModel to allocate
+ *
+ */
 void msvmmaj_allocate_model(struct MajModel *model)
 {
 	long n = model->n;
@@ -360,6 +517,16 @@ void msvmmaj_allocate_model(struct MajModel *model)
 
 }	
 
+/**
+ * @brief Free allocated MajModel struct
+ *
+ * @details
+ * Simply free a previously allocated MajModel by freeing all its component
+ * arrays. Note that the model struct itself is also freed here.
+ *
+ * @param[in] 	model 	MajModel to free
+ *
+ */
 void msvmmaj_free_model(struct MajModel *model)
 {
 	free(model->W);
@@ -376,10 +543,19 @@ void msvmmaj_free_model(struct MajModel *model)
 	free(model);
 }
 
+/**
+ * @brief Free allocated MajData struct 
+ *
+ * @details
+ * Simply free a previously allocated MajData struct by freeing all its
+ * components. Note that the data struct itself is also freed here.
+ *
+ * @param[in] 	data 	MajData struct to free
+ *
+ */
 void msvmmaj_free_data(struct MajData *data)
 {
 	free(data->Z);
 	free(data->y);
 	free(data);
 }
-