SHOGUN
v2.0.0
|
00001 /* 00002 * This program is free software; you can redistribute it and/or modify 00003 * it under the terms of the GNU General Public License as published by 00004 * the Free Software Foundation; either version 3 of the License, or 00005 * (at your option) any later version. 00006 * 00007 * Written (W) 2012 Viktor Gal 00008 * Copyright (C) 2007-12 Andrea Vedaldi and Brian Fulkerson. 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 /* precompute the sampled periodicized spectrum */ 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 /* fill table */ 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 } /* next mantissa */ 00159 } /* next exponent */ 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 /* throw exception */ 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 /* project on the postivie orthant (see PAMI) */ 00289 kappa_hat = CMath::max (kappa_hat, 0.0); 00290 break; 00291 default: 00292 /* throw exception */ 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 /* assert for in vector */ 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 /* break value into exponent and mantissa */ 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 //uint64_t featureDimension = 2*m_order+1; 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 /* register variables */ 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 }