Merge pull request #374 from vigsterkr/master

Add Homogeneous Kernel Mapping preprocessor
shogun-toolbox · Feb 17, 2012 · 9abe870 · 9abe870
2 parents d7d10b9 + 53cab65
commit 9abe870
Show file tree

Hide file tree

Showing 8 changed files with 474 additions and 5 deletions.
diff --git a/src/interfaces/modular/Preprocessor.i b/src/interfaces/modular/Preprocessor.i
@@ -16,6 +16,7 @@
 %rename(LogPlusOne) CLogPlusOne;
 %rename(PruneVarSubMean) CPruneVarSubMean;
 %rename(RandomFourierGaussPreproc) CRandomFourierGaussPreproc;
+%rename(HomogeneousKernelMap) CHomogeneousKernelMap;
 
 %rename(DimensionReductionPreprocessor) CDimensionReductionPreprocessor;
 %rename(PCA) CPCA;
@@ -98,6 +99,7 @@ namespace shogun
 %include <shogun/preprocessor/LogPlusOne.h>
 %include <shogun/preprocessor/PruneVarSubMean.h>
 %include <shogun/preprocessor/RandomFourierGaussPreproc.h>
+%include <shogun/preprocessor/HomogeneousKernelMap.h>
 
 %include <shogun/preprocessor/PCA.h>
 %include <shogun/preprocessor/KernelPCA.h>

diff --git a/src/interfaces/modular/Preprocessor_includes.i b/src/interfaces/modular/Preprocessor_includes.i
@@ -10,6 +10,7 @@
 #include <shogun/preprocessor/LogPlusOne.h>
 #include <shogun/preprocessor/PruneVarSubMean.h>
 #include <shogun/preprocessor/RandomFourierGaussPreproc.h>
+#include <shogun/preprocessor/HomogeneousKernelMap.h>
 
 #include <shogun/preprocessor/DimensionReductionPreprocessor.h>
 #include <shogun/preprocessor/PCA.h>

diff --git a/src/shogun/kernel/HistogramIntersectionKernel.cpp b/src/shogun/kernel/HistogramIntersectionKernel.cpp
@@ -31,7 +31,7 @@ CHistogramIntersectionKernel::CHistogramIntersectionKernel(int32_t size)
 CHistogramIntersectionKernel::CHistogramIntersectionKernel(
 	CSimpleFeatures<float64_t>* l, CSimpleFeatures<float64_t>* r,
 	float64_t beta, int32_t size)
-: CDotKernel(size), m_beta(1.0)
+: CDotKernel(size), m_beta(beta)
 {
 	init(l,r);
 	register_params();

diff --git a/src/shogun/kernel/JensenShannonKernel.cpp b/src/shogun/kernel/JensenShannonKernel.cpp
@@ -63,11 +63,11 @@ float64_t CJensenShannonKernel::compute(int32_t idx_a, int32_t idx_b)
 		float64_t a_i = 0, b_i = 0;
 		float64_t ab = avec[i]+bvec[i];
 		if (avec[i] != 0)
-			a_i = avec[i]/2 * CMath::log2(ab/avec[i]);
+			a_i = avec[i] * CMath::log2(ab/avec[i]);
 		if (bvec[i] != 0)
-			b_i = bvec[i]/2 * CMath::log2(ab/bvec[i]);
+			b_i = bvec[i] * CMath::log2(ab/bvec[i]);
 
-		result += a_i + b_i;
+		result += 0.5*(a_i + b_i);
 	}
 
 	((CSimpleFeatures<float64_t>*) lhs)->free_feature_vector(avec, idx_a, afree);

diff --git a/src/shogun/mathematics/Math.h b/src/shogun/mathematics/Math.h
@@ -421,6 +421,11 @@ class CMath : public CSGObject
 			return ::sin(x);
 		}
 
+		static inline float64_t cos(float64_t x)
+		{
+			return ::cos(x);
+		}
+
 		static float64_t area_under_curve(float64_t* xy, int32_t len, bool reversed)
 		{
 			ASSERT(len>0 && xy);

diff --git a/src/shogun/preprocessor/HomogeneousKernelMap.cpp b/src/shogun/preprocessor/HomogeneousKernelMap.cpp
@@ -0,0 +1,336 @@
+/*
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * Written (W) 2012 Viktor Gal
+ * Copyright (C) 2007-12 Andrea Vedaldi and Brian Fulkerson.
+ */
+
+#include <shogun/io/SGIO.h>
+#include <shogun/preprocessor/HomogeneousKernelMap.h>
+
+#include <math.h>
+#include <iostream>
+
+using namespace shogun;
+
+CHomogeneousKernelMap::CHomogeneousKernelMap () 
+	: CSimplePreprocessor<float64_t> ()
+{
+
+}
+
+CHomogeneousKernelMap::CHomogeneousKernelMap 
+	(HomogeneousKernelType kernel, HomogeneousKernelMapWindowType wType,
+	 float64_t gamma, uint64_t order, float64_t period)
+	: CSimplePreprocessor<float64_t> (),
+		m_kernel (kernel),
+		m_window (wType),
+		m_gamma (gamma),
+		m_period (period),
+		m_order (order)
+
+{
+	init ();
+}
+
+CHomogeneousKernelMap::~CHomogeneousKernelMap() {
+	m_table.destroy_vector ();
+}
+
+bool CHomogeneousKernelMap::init (CFeatures* features) {
+	ASSERT(features->get_feature_class()==C_SIMPLE);
+	ASSERT(features->get_feature_type()==F_DREAL);
+
+	return true;
+}
+
+void CHomogeneousKernelMap::cleanup () {
+	m_table.destroy_vector ();
+}
+
+
+void CHomogeneousKernelMap::init () {
+	SG_DEBUG ("Initialising homogeneous kernel map...\n");
+	ASSERT (m_gamma > 0) ;
+
+  ASSERT (m_kernel == HomogeneousKernelIntersection ||
+          m_kernel == HomogeneousKernelChi2 ||
+          m_kernel == HomogeneousKernelJS);
+
+  ASSERT (m_window == HomogeneousKernelMapWindowUniform ||
+          m_window == HomogeneousKernelMapWindowRectangular);
+
+	if (m_period < 0) {
+	  switch (m_window) {
+	  case HomogeneousKernelMapWindowUniform:
+	    switch (m_kernel) {
+		    case HomogeneousKernelChi2:         
+					m_period = 5.86 * CMath::sqrt (static_cast<float64_t> (m_order))  + 3.65; 
+					break;
+		    case HomogeneousKernelJS:
+		      m_period = 6.64 * CMath::sqrt (static_cast<float64_t> (m_order))  + 7.24; 
+					break;
+		    case HomogeneousKernelIntersection: 
+					m_period = 2.38 * CMath::log (m_order + 0.8) + 5.6;
+					break;
+	    }
+	    break;
+	  case HomogeneousKernelMapWindowRectangular:
+	    switch (m_kernel) {
+		    case HomogeneousKernelChi2:
+					m_period = 8.80 * CMath::sqrt (m_order + 4.44) - 12.6; 
+					break;
+		    case HomogeneousKernelJS:
+					m_period = 9.63 * CMath::sqrt (m_order + 1.00) - 2.93;
+					break;
+		    case HomogeneousKernelIntersection: 
+					m_period = 2.00 * CMath::log (m_order + 0.99) + 3.52;
+					break;
+	    }
+	    break;
+	  }
+	  m_period = CMath::max (m_period, 1.0) ;
+	}
+
+  m_numSubdivisions = 8 + 8*m_order;
+  m_subdivision = 1.0 / m_numSubdivisions;
+  m_minExponent = -20;
+  m_maxExponent = 8;
+
+  int tableHeight = 2*m_order + 1 ;
+  int tableWidth = m_numSubdivisions * (m_maxExponent - m_minExponent + 1);
+	size_t numElements = (tableHeight * tableWidth + 2*(1+m_order));
+	m_table = SGVector<float64_t> (NULL, numElements, true);
+	m_table.vlen = numElements;
+	m_table.vector = SG_CALLOC (float64_t, numElements);
+
+	/* check whether allocation was ok */
+	if (m_table.vector == NULL) {
+		throw ShogunException ("CHomogeneousKernelMap::init: could not allocate memory for vector");
+	}
+
+	int exponent;
+  uint64_t i = 0, j = 0;
+  float64_t* tablep = m_table.vector;
+  float64_t* kappa = m_table.vector + tableHeight * tableWidth;
+  float64_t* freq = kappa + (1+m_order);
+  float64_t L = 2.0 * CMath::PI / m_period;
+
+  /* precompute the sampled periodicized spectrum */
+  while (i <= m_order) {
+    freq[i] = j;
+    kappa[i] = get_smooth_spectrum (j * L);
+    ++ j;
+    if (kappa[i] > 0 || j >= 3*i) ++ i;
+  }
+
+  /* fill table */
+  for (exponent  = m_minExponent ;
+       exponent <= m_maxExponent ; ++ exponent) {
+
+    float64_t x, Lxgamma, Llogx, xgamma;
+    float64_t sqrt2kappaLxgamma;
+    float64_t mantissa = 1.0;
+
+    for (i = 0 ; i < m_numSubdivisions;
+         ++i, mantissa += m_subdivision) {
+      x = ldexp (mantissa, exponent);
+      xgamma = CMath::pow (x, m_gamma);
+      Lxgamma = L * xgamma;
+      Llogx = L * CMath::log (x);
+
+      *tablep++ = CMath::sqrt (Lxgamma * kappa[0]);
+      for (j = 1 ; j <= m_order; ++j) {
+        sqrt2kappaLxgamma = CMath::sqrt (2.0 * Lxgamma * kappa[j]);
+        *tablep++ = sqrt2kappaLxgamma * CMath::cos (freq[j] * Llogx);
+        *tablep++ = sqrt2kappaLxgamma * CMath::sin (freq[j] * Llogx);
+      }
+    } /* next mantissa */
+  } /* next exponent */
+
+	/* register variables */
+	m_parameters->add ((machine_int_t*) &m_kernel, "kernel", "Kernel type to use.");
+	m_parameters->add ((machine_int_t*) &m_window, "window", "Window type to use.");
+	m_parameters->add (&m_gamma, "gamma", "Homogeneity order.");
+	m_parameters->add (&m_period, "period", "Approximation order");
+	m_parameters->add (&m_numSubdivisions, "numSubdivisions", "The number of sublevels");
+	m_parameters->add (&m_subdivision, "subdivision", "subdivision.");
+	m_parameters->add (&m_order, "order", "The order");
+	m_parameters->add (&m_minExponent, "minExponent", "Minimum exponent");
+	m_parameters->add (&m_maxExponent, "maxExponent", "Maximum exponent");
+	m_parameters->add (&m_table, "table", "Lookup-table");
+}
+
+SGMatrix<float64_t> CHomogeneousKernelMap::apply_to_feature_matrix (CFeatures* features) {
+	CSimpleFeatures<float64_t>* simple_features = (CSimpleFeatures<float64_t>*)features;
+	int32_t num_vectors = simple_features->get_num_vectors ();
+	int32_t num_features = simple_features->get_num_features ();
+	SGMatrix<float64_t> result (num_features*(2*m_order+1), num_vectors);
+
+	for (int i = 0; i < num_vectors; ++i) {
+		SGVector<float64_t> v = simple_features->get_feature_vector (i);
+		SGVector<float64_t> col (result.get_column_vector (i), result.num_rows);
+		apply_to_vector (v, col);
+	}
+
+	return result;
+}
+
+/// apply preproc on single feature vector
+SGVector<float64_t> CHomogeneousKernelMap::apply_to_feature_vector (SGVector<float64_t> vector) {
+	uint64_t featureDimension = 2*m_order+1;
+	uint64_t m_target_dim = vector.vlen * featureDimension;
+	SGVector<float64_t> result = SGVector<float64_t> (m_target_dim);
+
+	apply_to_vector (vector, result);
+
+	return result;
+}
+
+void CHomogeneousKernelMap::setKernelType (HomogeneousKernelType k) {
+	m_kernel = k;
+}
+
+HomogeneousKernelType CHomogeneousKernelMap::getKernelType () const {
+	return m_kernel;
+}
+
+void CHomogeneousKernelMap::setWindowType (HomogeneousKernelMapWindowType w) {
+	m_window = w;
+}
+
+HomogeneousKernelMapWindowType CHomogeneousKernelMap::getWindowType () const {
+	return m_window;
+}
+
+void CHomogeneousKernelMap::setGamma (float64_t g) {
+	m_gamma = g;
+}
+
+float64_t CHomogeneousKernelMap::getGamma (float64_t g) const {
+	return m_gamma;
+}
+
+void CHomogeneousKernelMap::setOrder (uint64_t o) {
+	m_order = o;
+}
+
+uint64_t CHomogeneousKernelMap::getOrder () const {
+	return m_order;
+}
+
+void CHomogeneousKernelMap::setPeriod (float64_t p) {
+	m_period = p;
+}
+
+float64_t CHomogeneousKernelMap::getPeriod () const {
+	return m_period;
+}
+
+inline float64_t
+CHomogeneousKernelMap::get_spectrum (float64_t omega) const {
+  switch (m_kernel) {
+    case HomogeneousKernelIntersection:
+      return (2.0 / CMath::PI) / (1 + 4 * omega*omega);
+    case HomogeneousKernelChi2:
+      return 2.0 / (CMath::exp (CMath::PI * omega) + CMath::exp (-CMath::PI * omega)) ;
+    case HomogeneousKernelJS:
+      return (2.0 / CMath::log (4.0)) *
+      2.0 / (CMath::exp (CMath::PI * omega) + CMath::exp (-CMath::PI * omega)) /
+      (1 + 4 * omega*omega);
+    default:
+			/* throw exception */
+			throw ShogunException ("CHomogeneousKernelMap::get_spectrum: no valid kernel has been set!");
+  }
+}
+
+inline float64_t 
+CHomogeneousKernelMap::sinc (float64_t x) const {
+  if (x == 0.0) return 1.0 ;
+  return CMath::sin (x) / x;
+}
+
+inline float64_t
+CHomogeneousKernelMap::get_smooth_spectrum (float64_t omega) const {
+  float64_t kappa_hat = 0;
+  float64_t omegap ;
+  float64_t epsilon = 1e-2;
+  float64_t const omegaRange = 2.0 / (m_period * epsilon);
+  float64_t const domega = 2 * omegaRange / (2 * 1024.0 + 1);
+  switch (m_window) {
+    case HomogeneousKernelMapWindowUniform:
+      kappa_hat = get_spectrum (omega);
+      break;
+    case HomogeneousKernelMapWindowRectangular:
+      for (omegap = - omegaRange ; omegap <= omegaRange ; omegap += domega) {
+        float64_t win = sinc ((m_period/2.0) * omegap);
+        win *= (m_period/(2.0*CMath::PI));
+        kappa_hat += win * get_spectrum (omegap + omega);
+      }
+      kappa_hat *= domega;
+      /* project on the postivie orthant (see PAMI) */
+      kappa_hat = CMath::max (kappa_hat, 0.0);
+      break;
+    default:
+			/* throw exception */
+			throw ShogunException ("CHomogeneousKernelMap::get_smooth_spectrum: no valid kernel has been set!");
+  }
+  return kappa_hat;
+}
+
+inline void CHomogeneousKernelMap::apply_to_vector (const SGVector<float64_t>& in_v, 
+																										SGVector<float64_t>& out_v) const 
+{
+	/* assert for in and out vectors */
+	ASSERT (in_v.vlen > 0 && out_v.vlen);
+	ASSERT (in_v.vector != NULL && out_v.vector != NULL);
+
+	uint64_t featureDimension = 2*m_order+1;
+
+	for (int k = 0; k < in_v.vlen; ++k) {
+	  /* break value into exponent and mantissa */
+	  int exponent;
+	  int unsigned j;
+	  float64_t mantissa = frexp (in_v[k], &exponent);
+	  float64_t sign = (mantissa >= 0.0) ? +1.0 : -1.0;
+	  mantissa *= 2*sign;
+	  exponent -- ;
+
+	  if (mantissa == 0 ||
+	      exponent <= m_minExponent ||
+	      exponent >= m_maxExponent) 
+		{
+	    for (j = 0 ; j <= m_order ; ++j) {
+	      out_v[k*featureDimension+j] = 0.0;
+//				*destination = (T) 0.0 ;
+//	      destination += stride ;
+	    }
+	    continue;
+	  }
+
+	  //uint64_t featureDimension = 2*m_order+1;
+    float64_t const * v1 = m_table.vector +
+    	(exponent - m_minExponent) * m_numSubdivisions * featureDimension;
+    float64_t const * v2;
+    float64_t f1, f2;
+
+    mantissa -= 1.0;
+    while (mantissa >= m_subdivision) {
+      mantissa -= m_subdivision;
+      v1 += featureDimension;
+    }
+
+    v2 = v1 + featureDimension;
+    for (j = 0 ; j < featureDimension ; ++j) {
+      f1 = *v1++;
+      f2 = *v2++;
+
+			out_v[k*featureDimension+j] = sign * ((f2 - f1) * (m_numSubdivisions * mantissa) + f1);
+//      *destination = sign * ((f2 - f1) * (m_numSubdivisions * mantissa) + f1) ;
+//      destination += stride ;
+    }
+	}	
+}