1 files changed, 155 insertions, 39 deletions

diff --git a/src/msvmmaj_kernel.c b/src/msvmmaj_kernel.c
index 5ac138c..a2c1193 100644
--- a/src/msvmmaj_kernel.c
+++ b/src/msvmmaj_kernel.c
@@ -32,14 +32,60 @@
 void msvmmaj_make_kernel(struct MajModel *model, struct MajData *data)
 {
 	long i, j;
+	// Determine if a kernel needs to be computed. This is not the case if 
+	// a LINEAR kernel is requested in the model, or if the requested 
+	// kernel is already in the data. 
+
 	if (model->kerneltype == K_LINEAR) {
-		model->J = Calloc(double, model->m+1);
-		for (i=1; i<model->m+1; i++) 
-			matrix_set(model->J, 1, i, 0, 1.0);
+		data->J = Calloc(double, data->m+1);
+		for (i=1; i<data->m+1; i++) {
+			matrix_set(data->J, 1, i, 0, 1.0);
+		}
 		return;
 	}
 
-	long n = model->n;
+	/*
+	switch (model->kerneltype) {
+		case K_LINEAR:
+			// if data has another kernel, free that matrix and 
+			// assign Z to RAW
+			if (data->kerneltype != K_LINEAR) {
+				free(data->Z);
+				data->Z = data->RAW;
+			}
+			data->J = Calloc(double, data->m+1);
+			for (i=1; i<model->m+1; i++) {
+				matrix_set(data->J, 1, i, 0, 1.0);
+			}
+			return;
+		case K_POLY:
+			// if data has another kernel, we need to recalculate
+			if (data->kerneltype != K_POLY) {
+				break;
+			}
+			// if it is poly, we only recalculate if the kernel 
+			// parameters differ
+			if (data->kernelparam[0] == model->kernelparam[0] &&
+				data->kernelparam[1] == model->kernelparam[1] &&
+				data->kernelparam[2] == model->kernelparam[2])
+				// < do something with J ?
+				return;
+		case K_RBF:
+			if (data->kerneltype != K_RBF)
+				break;
+			if (data->kernelparam[0] == model->kernelparam[0])
+				// < do something with J ?
+				return;
+		case K_SIGMOID:
+			if (data->kerneltype != K_SIGMOID)
+				break;
+			if (data->kernelparam[0] == model->kernelparam[0] &&
+				data->kernelparam[1] == model->kernelparam[1])
+				// < do something with J ?
+				return;
+	}
+	*/
+	long n = data->n;
 	double value;
 	double *x1, *x2;
 	double *K = Calloc(double, n*n);
@@ -55,7 +101,7 @@ void msvmmaj_make_kernel(struct MajModel *model, struct MajData *data)
 				value = msvmmaj_compute_rbf(x1, x2, 
 						model->kernelparam, data->m);
 			else if (model->kerneltype == K_SIGMOID)
-				value = msvmmaj_compute_rbf(x1, x2, 
+				value = msvmmaj_compute_sigmoid(x1, x2, 
 						model->kernelparam, data->m);
 			else {
 				fprintf(stderr, "Unknown kernel type in "
@@ -67,25 +113,28 @@ void msvmmaj_make_kernel(struct MajModel *model, struct MajData *data)
 		}
 	}
 
-	double *P = Malloc(double, n*n);
-	double *Lambda = Malloc(double, n);
-	long num_eigen = msvmmaj_make_eigen(K, n, P, Lambda);
+	double *P = NULL;
+	double *Sigma = NULL;
+	long num_eigen = msvmmaj_make_eigen(K, n, &P, &Sigma);
+	//printf("num eigen: %li\n", num_eigen);
+	data->m = num_eigen;
+	model->m = num_eigen;
 
 	// copy eigendecomp to data
-	data->Z = realloc(data->Z, n*(n+1)*sizeof(double));
+	data->Z = Calloc(double, n*(num_eigen+1));
 	for (i=0; i<n; i++) {
-		for (j=0; j<n; j++)
-			matrix_set(data->Z, n+1, i, j+1,
-					matrix_get(P, n, i, j));
-		matrix_set(data->Z, n+1, i, 0, 1.0);
+		for (j=0; j<num_eigen; j++) {
+			value = matrix_get(P, num_eigen, i, j);
+			matrix_set(data->Z, num_eigen+1, i, j, value);
+		}
+		matrix_set(data->Z, num_eigen+1, i, 0, 1.0);
 	}
-	data->m = n;
 
 	// Set the regularization matrix (change if not full rank used)
-	model->J = Calloc(double, model->m+1);
-	for (i=1; i<model->m+1; i++) {
-		value = 1.0/matrix_get(Lambda, 1, i-1, 0);
-		matrix_set(model->J, 1, i, 0, value);
+	data->J = Calloc(double, data->m+1);
+	for (i=1; i<data->m+1; i++) {
+		value = 1.0/matrix_get(Sigma, 1, i-1, 0);
+		matrix_set(data->J, 1, i, 0, value);
 	}
 
 	// let data know what it's made of
@@ -109,23 +158,25 @@ void msvmmaj_make_kernel(struct MajModel *model, struct MajData *data)
 			data->kernelparam[0] = model->kernelparam[0];
 			data->kernelparam[1] = model->kernelparam[1];
 	}
-	model->m = n;
 	free(K);
+	free(Sigma);
+	free(P);
 }
 
 /**
- * @brief Find the eigendecomposition of a kernel matrix.
+ * @brief Find the (reduced) eigendecomposition of a kernel matrix.
  *
  * @details.
  * tbd
  *
- *
  */
-long msvmmaj_make_eigen(double *K, long n, double *P, double *Lambda)
+long msvmmaj_make_eigen(double *K, long n, double **P, double **Sigma)
 {
 	int M, status, LWORK, *IWORK, *IFAIL;
-	double abstol, *WORK;
+	long i, j, num_eigen, cutoff_idx;
+	double max_eigen, abstol, *WORK;
 
+	double *tempSigma = Malloc(double, n);
 	double *tempP = Malloc(double, n*n);
 
 	IWORK = Malloc(int, 5*n);
@@ -150,7 +201,7 @@ long msvmmaj_make_eigen(double *K, long n, double *P, double *Lambda)
 			0,
 			abstol,
 			&M,
-			Lambda,
+			tempSigma,
 			tempP,
 			n,
 			WORK,
@@ -174,7 +225,7 @@ long msvmmaj_make_eigen(double *K, long n, double *P, double *Lambda)
 			0,
 			abstol,
 			&M,
-			Lambda,
+			tempSigma,
 			tempP,
 			n,
 			WORK,
@@ -182,29 +233,94 @@ long msvmmaj_make_eigen(double *K, long n, double *P, double *Lambda)
 			IWORK,
 			IFAIL);
 
-	printf("status = %i\n", status);
-	printf("Number of eigenvalues found: %i\n", M);
-
 	if (status != 0) {
 		fprintf(stderr, "Nonzero exit status from dsyevx. Exiting...");
 		exit(1);
 	}
-	// Here you can put code to select the necessary eigenvectors, 
-	// depending on the size of the eigenvalues.
-	// For now, let's just print the eigenvalues and exit
+
+	// Select the desired number of eigenvalues, depending on their size.  
+	// dsyevx sorts eigenvalues in ascending order.
+	//
+	max_eigen = tempSigma[n-1];
+	cutoff_idx = 0;
+
+	for (i=0; i<n; i++)
+		if (tempSigma[i]/max_eigen > 1e-10 ) {
+			cutoff_idx = i;
+			break;
+		}
+
+	num_eigen = n - cutoff_idx;
 	
-	print_matrix(Lambda, n, 1);
+	*Sigma = Calloc(double, num_eigen);
 	
-	// revert P to row-major order	
-	long i, j;
-	for (i=0; i<n; i++)
-		for (j=0; j<n; j++)
-			P[i*n+j] = tempP[j*n+i];
+	for (i=0; i<num_eigen; i++) {
+		(*Sigma)[i] = tempSigma[n-1 - i];
+	}
 
+	// revert P to row-major order and copy only the the columns 
+	// corresponding to the selected eigenvalues
+	//
+	*P = Calloc(double, n*num_eigen);	
+	for (j=n-1; j>n-1-num_eigen; j--) {
+		for (i=0; i<n; i++) {
+			(*P)[i*num_eigen + (n-1)-j] = tempP[i + j*n];
+		}
+	}
+
+	free(tempSigma);	
 	free(tempP);
 
-	// replace by number of columns of P
-	return n;
+	return num_eigen;
+}
+
+void msvmmaj_make_crosskernel(struct MajModel *model,
+	       	struct MajData *data_train, struct MajData *data_test,
+	       	double **K2)
+{
+	long i, j;
+	long n_train = data_train->n;
+	long n_test = data_test->n;
+	long m = data_test->m;
+	double value;
+	double *x1, *x2;
+
+	*K2 = Calloc(double, n_test*n_train);
+
+	//printf("Training RAW\n");
+	//print_matrix(data_train->RAW, n_train, m+1);
+
+	//printf("Testing RAW\n");
+	//print_matrix(data_test->RAW, n_test, m+1);
+
+	for (i=0; i<n_test; i++) {
+		for (j=0; j<n_train; j++) {
+			x1 = &data_test->RAW[i*(m+1)+1];
+			x2 = &data_train->RAW[j*(m+1)+1];
+			if (model->kerneltype == K_POLY)
+				value = msvmmaj_compute_poly(x1, x2,
+						model->kernelparam,
+					       	m);
+			else if (model->kerneltype == K_RBF)
+				value = msvmmaj_compute_rbf(x1, x2,
+						model->kernelparam,
+					       	m);
+			else if (model->kerneltype == K_SIGMOID)
+				value = msvmmaj_compute_sigmoid(x1, x2,
+						model->kernelparam,
+					       	m);
+			else {
+				fprintf(stderr, "Unknown kernel type in "
+						"msvmmaj_make_crosskernel\n");
+				exit(1);
+			}
+			matrix_set((*K2), n_train, i, j, value);
+		}
+	}
+
+	//printf("cross K2:\n");
+	//print_matrix((*K2), n_test, n_train);
+
 }
 
 /**