SHOGUN  3.2.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
HomogeneousKernelMap.cpp
Go to the documentation of this file.
1 /*
2  * This program is free software; you can redistribute it and/or modify
3  * it under the terms of the GNU General Public License as published by
4  * the Free Software Foundation; either version 3 of the License, or
5  * (at your option) any later version.
6  *
7  * Written (W) 2012 Viktor Gal
8  * Copyright (C) 2007-12 Andrea Vedaldi and Brian Fulkerson.
9  */
10 
11 #include <shogun/io/SGIO.h>
13 
14 #include <math.h>
15 #include <iostream>
16 
17 using namespace shogun;
18 
23  m_gamma(1.0),
24  m_period(-1),
25  m_order(1)
26 {
27  init ();
28  register_params ();
29 }
30 
33  uint64_t order, float64_t period)
35  m_kernel(kernel),
36  m_window(wType),
37  m_gamma(gamma),
38  m_period(period),
39  m_order(order)
40 
41 {
42  init ();
43  register_params ();
44 }
45 
47 {
48 }
49 
50 bool CHomogeneousKernelMap::init(CFeatures* features)
51 {
52  ASSERT(features->get_feature_class()==C_DENSE)
53  ASSERT(features->get_feature_type()==F_DREAL)
54 
55  return true;
56 }
57 
59 {
60  m_table=SGVector<float64_t>();
61 }
62 
63 
64 void CHomogeneousKernelMap::init()
65 {
66  SG_DEBUG ("Initialising homogeneous kernel map...\n")
67  ASSERT (m_gamma > 0)
68 
70  m_kernel == HomogeneousKernelChi2 ||
71  m_kernel == HomogeneousKernelJS);
72 
75 
76  if (m_period < 0) {
77  switch (m_window) {
79  switch (m_kernel) {
81  m_period = 5.86 * CMath::sqrt (static_cast<float64_t> (m_order)) + 3.65;
82  break;
83  case HomogeneousKernelJS:
84  m_period = 6.64 * CMath::sqrt (static_cast<float64_t> (m_order)) + 7.24;
85  break;
87  m_period = 2.38 * CMath::log (m_order + 0.8) + 5.6;
88  break;
89  }
90  break;
91  case HomogeneousKernelMapWindowRectangular:
92  switch (m_kernel) {
94  m_period = 8.80 * CMath::sqrt (m_order + 4.44) - 12.6;
95  break;
96  case HomogeneousKernelJS:
97  m_period = 9.63 * CMath::sqrt (m_order + 1.00) - 2.93;
98  break;
100  m_period = 2.00 * CMath::log (m_order + 0.99) + 3.52;
101  break;
102  }
103  break;
104  }
105  m_period = CMath::max (m_period, 1.0) ;
106  }
107 
108  m_numSubdivisions = 8 + 8*m_order;
109  m_subdivision = 1.0 / m_numSubdivisions;
110  m_minExponent = -20;
111  m_maxExponent = 8;
112 
113  int tableHeight = 2*m_order + 1 ;
114  int tableWidth = m_numSubdivisions * (m_maxExponent - m_minExponent + 1);
115  size_t numElements = (tableHeight * tableWidth + 2*(1+m_order));
116  if (unsigned(m_table.vlen) != numElements) {
117  SG_DEBUG ("reallocating... %d -> %d\n", m_table.vlen, numElements)
118  m_table.vector = SG_REALLOC (float64_t, m_table.vector, m_table.vlen, numElements);
119  m_table.vlen = numElements;
120  }
121 
122  int exponent;
123  uint64_t i = 0, j = 0;
124  float64_t* tablep = m_table.vector;
125  float64_t* kappa = m_table.vector + tableHeight * tableWidth;
126  float64_t* freq = kappa + (1+m_order);
127  float64_t L = 2.0 * CMath::PI / m_period;
128 
129  /* precompute the sampled periodicized spectrum */
130  while (i <= m_order) {
131  freq[i] = j;
132  kappa[i] = get_smooth_spectrum (j * L);
133  ++ j;
134  if (kappa[i] > 0 || j >= 3*i) ++ i;
135  }
136 
137  /* fill table */
138  for (exponent = m_minExponent ;
139  exponent <= m_maxExponent ; ++ exponent) {
140 
141  float64_t x, Lxgamma, Llogx, xgamma;
142  float64_t sqrt2kappaLxgamma;
143  float64_t mantissa = 1.0;
144 
145  for (i = 0 ; i < m_numSubdivisions;
146  ++i, mantissa += m_subdivision) {
147  x = ldexp (mantissa, exponent);
148  xgamma = CMath::pow (x, m_gamma);
149  Lxgamma = L * xgamma;
150  Llogx = L * CMath::log (x);
151 
152  *tablep++ = CMath::sqrt (Lxgamma * kappa[0]);
153  for (j = 1 ; j <= m_order; ++j) {
154  sqrt2kappaLxgamma = CMath::sqrt (2.0 * Lxgamma * kappa[j]);
155  *tablep++ = sqrt2kappaLxgamma * CMath::cos (freq[j] * Llogx);
156  *tablep++ = sqrt2kappaLxgamma * CMath::sin (freq[j] * Llogx);
157  }
158  } /* next mantissa */
159  } /* next exponent */
160 
161 }
162 
164 {
165  CDenseFeatures<float64_t>* simple_features = (CDenseFeatures<float64_t>*)features;
166  int32_t num_vectors = simple_features->get_num_vectors ();
167  int32_t num_features = simple_features->get_num_features ();
168 
169  SGMatrix<float64_t> feature_matrix(num_features*(2*m_order+1),num_vectors);
170  for (int i = 0; i < num_vectors; ++i)
171  {
172  SGVector<float64_t> transformed = apply_to_vector(simple_features->get_feature_vector(i));
173  for (int j=0; j<transformed.vlen; j++)
174  feature_matrix(j,i) = transformed[j];
175  }
176 
177  simple_features->set_feature_matrix(feature_matrix);
178 
179  return feature_matrix;
180 }
181 
184 {
185  SGVector<float64_t> result = apply_to_vector(vector);
186  return result;
187 }
188 
190 {
191  m_kernel = k;
192  init ();
193 }
194 
196 {
197  return m_kernel;
198 }
199 
201 {
202  m_window = w;
203  init ();
204 }
205 
207 {
208  return m_window;
209 }
210 
212 {
213  m_gamma = g;
214  init ();
215 }
216 
217 float64_t CHomogeneousKernelMap::get_gamma(float64_t g) const
218 {
219  return m_gamma;
220 }
221 
223 {
224  m_order = o;
225  init ();
226 }
227 
229 {
230  return m_order;
231 }
232 
234 {
235  m_period = p;
236  init ();
237 }
238 
240 {
241  return m_period;
242 }
243 
244 inline float64_t
245 CHomogeneousKernelMap::get_spectrum(float64_t omega) const
246 {
247  switch (m_kernel) {
249  return (2.0 / CMath::PI) / (1 + 4 * omega*omega);
251  return 2.0 / (CMath::exp (CMath::PI * omega) + CMath::exp (-CMath::PI * omega)) ;
252  case HomogeneousKernelJS:
253  return (2.0 / CMath::log (4.0)) *
254  2.0 / (CMath::exp (CMath::PI * omega) + CMath::exp (-CMath::PI * omega)) /
255  (1 + 4 * omega*omega);
256  default:
257  /* throw exception */
258  throw ShogunException ("CHomogeneousKernelMap::get_spectrum: no valid kernel has been set!");
259  }
260 }
261 
262 inline float64_t
263 CHomogeneousKernelMap::sinc(float64_t x) const
264 {
265  if (x == 0.0) return 1.0 ;
266  return CMath::sin (x) / x;
267 }
268 
269 inline float64_t
270 CHomogeneousKernelMap::get_smooth_spectrum(float64_t omega) const
271 {
272  float64_t kappa_hat = 0;
273  float64_t omegap ;
274  float64_t epsilon = 1e-2;
275  float64_t const omegaRange = 2.0 / (m_period * epsilon);
276  float64_t const domega = 2 * omegaRange / (2 * 1024.0 + 1);
277  switch (m_window) {
279  kappa_hat = get_spectrum (omega);
280  break;
281  case HomogeneousKernelMapWindowRectangular:
282  for (omegap = - omegaRange ; omegap <= omegaRange ; omegap += domega) {
283  float64_t win = sinc ((m_period/2.0) * omegap);
284  win *= (m_period/(2.0*CMath::PI));
285  kappa_hat += win * get_spectrum (omegap + omega);
286  }
287  kappa_hat *= domega;
288  /* project on the postivie orthant (see PAMI) */
289  kappa_hat = CMath::max (kappa_hat, 0.0);
290  break;
291  default:
292  /* throw exception */
293  throw ShogunException ("CHomogeneousKernelMap::get_smooth_spectrum: no valid kernel has been set!");
294  }
295  return kappa_hat;
296 }
297 
298 SGVector<float64_t> CHomogeneousKernelMap::apply_to_vector(const SGVector<float64_t>& in_v) const
299 {
300  /* assert for in vector */
301  ASSERT (in_v.vlen > 0)
302  ASSERT (in_v.vector != NULL)
303 
304  uint64_t featureDimension = 2*m_order+1;
305 
306  SGVector<float64_t> out_v(featureDimension*in_v.vlen);
307 
308  for (int k = 0; k < in_v.vlen; ++k) {
309  /* break value into exponent and mantissa */
310  int exponent;
311  int unsigned j;
312  float64_t mantissa = frexp (in_v[k], &exponent);
313  float64_t sign = (mantissa >= 0.0) ? +1.0 : -1.0;
314  mantissa *= 2*sign;
315  exponent -- ;
316 
317  if (mantissa == 0 ||
318  exponent <= m_minExponent ||
319  exponent >= m_maxExponent)
320  {
321  for (j = 0 ; j <= m_order ; ++j) {
322  out_v[k*featureDimension+j] = 0.0;
323  }
324  continue;
325  }
326 
327  //uint64_t featureDimension = 2*m_order+1;
328  float64_t const * v1 = m_table.vector +
329  (exponent - m_minExponent) * m_numSubdivisions * featureDimension;
330  float64_t const * v2;
331  float64_t f1, f2;
332 
333  mantissa -= 1.0;
334  while (mantissa >= m_subdivision) {
335  mantissa -= m_subdivision;
336  v1 += featureDimension;
337  }
338 
339  v2 = v1 + featureDimension;
340  for (j = 0 ; j < featureDimension ; ++j) {
341  f1 = *v1++;
342  f2 = *v2++;
343 
344  out_v[k*featureDimension+j] = sign * ((f2 - f1) * (m_numSubdivisions * mantissa) + f1);
345  }
346  }
347  return out_v;
348 }
349 
350 void CHomogeneousKernelMap::register_params()
351 {
352  /* register variables */
353  SG_ADD((machine_int_t*) &m_kernel, "kernel", "Kernel type to use.",MS_AVAILABLE);
354  SG_ADD((machine_int_t*) &m_window, "window", "Window type to use.",MS_AVAILABLE);
355  SG_ADD(&m_gamma, "gamma", "Homogeneity order.",MS_AVAILABLE);
356  SG_ADD(&m_period, "period", "Approximation order",MS_NOT_AVAILABLE);
357  SG_ADD(&m_numSubdivisions, "num_subdivisions", "The number of sublevels",MS_NOT_AVAILABLE);
358  SG_ADD(&m_subdivision, "subdivision", "subdivision.",MS_NOT_AVAILABLE);
359  SG_ADD(&m_order, "order", "The order",MS_AVAILABLE);
360  SG_ADD(&m_minExponent, "min_exponent", "Minimum exponent",MS_NOT_AVAILABLE);
361  SG_ADD(&m_maxExponent, "max_exponent", "Maximum exponent",MS_NOT_AVAILABLE);
362  SG_ADD(&m_table, "table", "Lookup-table",MS_NOT_AVAILABLE);
363 }

SHOGUN Machine Learning Toolbox - Documentation