00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #include <shogun/io/SGIO.h>
00012 #include <shogun/preprocessor/HomogeneousKernelMap.h>
00013
00014 #include <math.h>
00015 #include <iostream>
00016
00017 using namespace shogun;
00018
00019 CHomogeneousKernelMap::CHomogeneousKernelMap()
00020 : CDensePreprocessor<float64_t>(),
00021 m_kernel(HomogeneousKernelIntersection),
00022 m_window(HomogeneousKernelMapWindowRectangular),
00023 m_gamma(1.0),
00024 m_period(-1),
00025 m_order(1)
00026 {
00027 init ();
00028 register_params ();
00029 }
00030
00031 CHomogeneousKernelMap::CHomogeneousKernelMap(HomogeneousKernelType kernel,
00032 HomogeneousKernelMapWindowType wType, float64_t gamma,
00033 uint64_t order, float64_t period)
00034 : CDensePreprocessor<float64_t>(),
00035 m_kernel(kernel),
00036 m_window(wType),
00037 m_gamma(gamma),
00038 m_period(period),
00039 m_order(order)
00040
00041 {
00042 init ();
00043 register_params ();
00044 }
00045
00046 CHomogeneousKernelMap::~CHomogeneousKernelMap()
00047 {
00048 }
00049
00050 bool CHomogeneousKernelMap::init(CFeatures* features)
00051 {
00052 ASSERT(features->get_feature_class()==C_DENSE);
00053 ASSERT(features->get_feature_type()==F_DREAL);
00054
00055 return true;
00056 }
00057
00058 void CHomogeneousKernelMap::cleanup()
00059 {
00060 m_table=SGVector<float64_t>();
00061 }
00062
00063
00064 void CHomogeneousKernelMap::init()
00065 {
00066 SG_DEBUG ("Initialising homogeneous kernel map...\n");
00067 ASSERT (m_gamma > 0) ;
00068
00069 ASSERT (m_kernel == HomogeneousKernelIntersection ||
00070 m_kernel == HomogeneousKernelChi2 ||
00071 m_kernel == HomogeneousKernelJS);
00072
00073 ASSERT (m_window == HomogeneousKernelMapWindowUniform ||
00074 m_window == HomogeneousKernelMapWindowRectangular);
00075
00076 if (m_period < 0) {
00077 switch (m_window) {
00078 case HomogeneousKernelMapWindowUniform:
00079 switch (m_kernel) {
00080 case HomogeneousKernelChi2:
00081 m_period = 5.86 * CMath::sqrt (static_cast<float64_t> (m_order)) + 3.65;
00082 break;
00083 case HomogeneousKernelJS:
00084 m_period = 6.64 * CMath::sqrt (static_cast<float64_t> (m_order)) + 7.24;
00085 break;
00086 case HomogeneousKernelIntersection:
00087 m_period = 2.38 * CMath::log (m_order + 0.8) + 5.6;
00088 break;
00089 }
00090 break;
00091 case HomogeneousKernelMapWindowRectangular:
00092 switch (m_kernel) {
00093 case HomogeneousKernelChi2:
00094 m_period = 8.80 * CMath::sqrt (m_order + 4.44) - 12.6;
00095 break;
00096 case HomogeneousKernelJS:
00097 m_period = 9.63 * CMath::sqrt (m_order + 1.00) - 2.93;
00098 break;
00099 case HomogeneousKernelIntersection:
00100 m_period = 2.00 * CMath::log (m_order + 0.99) + 3.52;
00101 break;
00102 }
00103 break;
00104 }
00105 m_period = CMath::max (m_period, 1.0) ;
00106 }
00107
00108 m_numSubdivisions = 8 + 8*m_order;
00109 m_subdivision = 1.0 / m_numSubdivisions;
00110 m_minExponent = -20;
00111 m_maxExponent = 8;
00112
00113 int tableHeight = 2*m_order + 1 ;
00114 int tableWidth = m_numSubdivisions * (m_maxExponent - m_minExponent + 1);
00115 size_t numElements = (tableHeight * tableWidth + 2*(1+m_order));
00116 if (unsigned(m_table.vlen) != numElements) {
00117 SG_DEBUG ("reallocating... %d -> %d\n", m_table.vlen, numElements);
00118 m_table.vector = SG_REALLOC (float64_t, m_table.vector, numElements);
00119 m_table.vlen = numElements;
00120 }
00121
00122 int exponent;
00123 uint64_t i = 0, j = 0;
00124 float64_t* tablep = m_table.vector;
00125 float64_t* kappa = m_table.vector + tableHeight * tableWidth;
00126 float64_t* freq = kappa + (1+m_order);
00127 float64_t L = 2.0 * CMath::PI / m_period;
00128
00129
00130 while (i <= m_order) {
00131 freq[i] = j;
00132 kappa[i] = get_smooth_spectrum (j * L);
00133 ++ j;
00134 if (kappa[i] > 0 || j >= 3*i) ++ i;
00135 }
00136
00137
00138 for (exponent = m_minExponent ;
00139 exponent <= m_maxExponent ; ++ exponent) {
00140
00141 float64_t x, Lxgamma, Llogx, xgamma;
00142 float64_t sqrt2kappaLxgamma;
00143 float64_t mantissa = 1.0;
00144
00145 for (i = 0 ; i < m_numSubdivisions;
00146 ++i, mantissa += m_subdivision) {
00147 x = ldexp (mantissa, exponent);
00148 xgamma = CMath::pow (x, m_gamma);
00149 Lxgamma = L * xgamma;
00150 Llogx = L * CMath::log (x);
00151
00152 *tablep++ = CMath::sqrt (Lxgamma * kappa[0]);
00153 for (j = 1 ; j <= m_order; ++j) {
00154 sqrt2kappaLxgamma = CMath::sqrt (2.0 * Lxgamma * kappa[j]);
00155 *tablep++ = sqrt2kappaLxgamma * CMath::cos (freq[j] * Llogx);
00156 *tablep++ = sqrt2kappaLxgamma * CMath::sin (freq[j] * Llogx);
00157 }
00158 }
00159 }
00160
00161 }
00162
00163 SGMatrix<float64_t> CHomogeneousKernelMap::apply_to_feature_matrix (CFeatures* features)
00164 {
00165 CDenseFeatures<float64_t>* simple_features = (CDenseFeatures<float64_t>*)features;
00166 int32_t num_vectors = simple_features->get_num_vectors ();
00167 int32_t num_features = simple_features->get_num_features ();
00168
00169 SGMatrix<float64_t> feature_matrix(num_features*(2*m_order+1),num_vectors);
00170 for (int i = 0; i < num_vectors; ++i)
00171 {
00172 SGVector<float64_t> transformed = apply_to_vector(simple_features->get_feature_vector(i));
00173 for (int j=0; j<transformed.vlen; j++)
00174 feature_matrix(j,i) = transformed[j];
00175 }
00176
00177 simple_features->set_feature_matrix(feature_matrix);
00178
00179 return feature_matrix;
00180 }
00181
00183 SGVector<float64_t> CHomogeneousKernelMap::apply_to_feature_vector(SGVector<float64_t> vector)
00184 {
00185 SGVector<float64_t> result = apply_to_vector(vector);
00186 return result;
00187 }
00188
00189 void CHomogeneousKernelMap::set_kernel_type(HomogeneousKernelType k)
00190 {
00191 m_kernel = k;
00192 init ();
00193 }
00194
00195 HomogeneousKernelType CHomogeneousKernelMap::get_kernel_type() const
00196 {
00197 return m_kernel;
00198 }
00199
00200 void CHomogeneousKernelMap::set_window_type(HomogeneousKernelMapWindowType w)
00201 {
00202 m_window = w;
00203 init ();
00204 }
00205
00206 HomogeneousKernelMapWindowType CHomogeneousKernelMap::get_window_type() const
00207 {
00208 return m_window;
00209 }
00210
00211 void CHomogeneousKernelMap::set_gamma(float64_t g)
00212 {
00213 m_gamma = g;
00214 init ();
00215 }
00216
00217 float64_t CHomogeneousKernelMap::get_gamma(float64_t g) const
00218 {
00219 return m_gamma;
00220 }
00221
00222 void CHomogeneousKernelMap::set_order(uint64_t o)
00223 {
00224 m_order = o;
00225 init ();
00226 }
00227
00228 uint64_t CHomogeneousKernelMap::get_order() const
00229 {
00230 return m_order;
00231 }
00232
00233 void CHomogeneousKernelMap::set_period(float64_t p)
00234 {
00235 m_period = p;
00236 init ();
00237 }
00238
00239 float64_t CHomogeneousKernelMap::get_period() const
00240 {
00241 return m_period;
00242 }
00243
00244 inline float64_t
00245 CHomogeneousKernelMap::get_spectrum(float64_t omega) const
00246 {
00247 switch (m_kernel) {
00248 case HomogeneousKernelIntersection:
00249 return (2.0 / CMath::PI) / (1 + 4 * omega*omega);
00250 case HomogeneousKernelChi2:
00251 return 2.0 / (CMath::exp (CMath::PI * omega) + CMath::exp (-CMath::PI * omega)) ;
00252 case HomogeneousKernelJS:
00253 return (2.0 / CMath::log (4.0)) *
00254 2.0 / (CMath::exp (CMath::PI * omega) + CMath::exp (-CMath::PI * omega)) /
00255 (1 + 4 * omega*omega);
00256 default:
00257
00258 throw ShogunException ("CHomogeneousKernelMap::get_spectrum: no valid kernel has been set!");
00259 }
00260 }
00261
00262 inline float64_t
00263 CHomogeneousKernelMap::sinc(float64_t x) const
00264 {
00265 if (x == 0.0) return 1.0 ;
00266 return CMath::sin (x) / x;
00267 }
00268
00269 inline float64_t
00270 CHomogeneousKernelMap::get_smooth_spectrum(float64_t omega) const
00271 {
00272 float64_t kappa_hat = 0;
00273 float64_t omegap ;
00274 float64_t epsilon = 1e-2;
00275 float64_t const omegaRange = 2.0 / (m_period * epsilon);
00276 float64_t const domega = 2 * omegaRange / (2 * 1024.0 + 1);
00277 switch (m_window) {
00278 case HomogeneousKernelMapWindowUniform:
00279 kappa_hat = get_spectrum (omega);
00280 break;
00281 case HomogeneousKernelMapWindowRectangular:
00282 for (omegap = - omegaRange ; omegap <= omegaRange ; omegap += domega) {
00283 float64_t win = sinc ((m_period/2.0) * omegap);
00284 win *= (m_period/(2.0*CMath::PI));
00285 kappa_hat += win * get_spectrum (omegap + omega);
00286 }
00287 kappa_hat *= domega;
00288
00289 kappa_hat = CMath::max (kappa_hat, 0.0);
00290 break;
00291 default:
00292
00293 throw ShogunException ("CHomogeneousKernelMap::get_smooth_spectrum: no valid kernel has been set!");
00294 }
00295 return kappa_hat;
00296 }
00297
00298 SGVector<float64_t> CHomogeneousKernelMap::apply_to_vector(const SGVector<float64_t>& in_v) const
00299 {
00300
00301 ASSERT (in_v.vlen > 0);
00302 ASSERT (in_v.vector != NULL);
00303
00304 uint64_t featureDimension = 2*m_order+1;
00305
00306 SGVector<float64_t> out_v(featureDimension*in_v.vlen);
00307
00308 for (int k = 0; k < in_v.vlen; ++k) {
00309
00310 int exponent;
00311 int unsigned j;
00312 float64_t mantissa = frexp (in_v[k], &exponent);
00313 float64_t sign = (mantissa >= 0.0) ? +1.0 : -1.0;
00314 mantissa *= 2*sign;
00315 exponent -- ;
00316
00317 if (mantissa == 0 ||
00318 exponent <= m_minExponent ||
00319 exponent >= m_maxExponent)
00320 {
00321 for (j = 0 ; j <= m_order ; ++j) {
00322 out_v[k*featureDimension+j] = 0.0;
00323 }
00324 continue;
00325 }
00326
00327
00328 float64_t const * v1 = m_table.vector +
00329 (exponent - m_minExponent) * m_numSubdivisions * featureDimension;
00330 float64_t const * v2;
00331 float64_t f1, f2;
00332
00333 mantissa -= 1.0;
00334 while (mantissa >= m_subdivision) {
00335 mantissa -= m_subdivision;
00336 v1 += featureDimension;
00337 }
00338
00339 v2 = v1 + featureDimension;
00340 for (j = 0 ; j < featureDimension ; ++j) {
00341 f1 = *v1++;
00342 f2 = *v2++;
00343
00344 out_v[k*featureDimension+j] = sign * ((f2 - f1) * (m_numSubdivisions * mantissa) + f1);
00345 }
00346 }
00347 return out_v;
00348 }
00349
00350 void CHomogeneousKernelMap::register_params()
00351 {
00352
00353 SG_ADD((machine_int_t*) &m_kernel, "kernel", "Kernel type to use.",MS_AVAILABLE);
00354 SG_ADD((machine_int_t*) &m_window, "window", "Window type to use.",MS_AVAILABLE);
00355 SG_ADD(&m_gamma, "gamma", "Homogeneity order.",MS_AVAILABLE);
00356 SG_ADD(&m_period, "period", "Approximation order",MS_NOT_AVAILABLE);
00357 SG_ADD(&m_numSubdivisions, "numSubdivisions", "The number of sublevels",MS_NOT_AVAILABLE);
00358 SG_ADD(&m_subdivision, "subdivision", "subdivision.",MS_NOT_AVAILABLE);
00359 SG_ADD(&m_order, "order", "The order",MS_AVAILABLE);
00360 SG_ADD(&m_minExponent, "minExponent", "Minimum exponent",MS_NOT_AVAILABLE);
00361 SG_ADD(&m_maxExponent, "maxExponent", "Maximum exponent",MS_NOT_AVAILABLE);
00362 SG_ADD(&m_table, "table", "Lookup-table",MS_NOT_AVAILABLE);
00363 }