SHOGUN  v1.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
DiffusionMaps.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) 2011 Sergey Lisitsyn
8  * Copyright (C) 2011 Berlin Institute of Technology and Max-Planck-Society
9  */
10 
13 #include <shogun/lib/config.h>
14 #ifdef HAVE_LAPACK
17 #include <shogun/io/SGIO.h>
18 #include <shogun/kernel/Kernel.h>
19 #include <shogun/lib/Signal.h>
20 #include <shogun/lib/Time.h>
21 
22 using namespace shogun;
23 
26 {
27  m_t = 10;
28 
29  init();
30 }
31 
33 {
34  m_parameters->add(&m_t, "t", "number of steps");
35 }
36 
38 {
39 }
40 
41 void CDiffusionMaps::set_t(int32_t t)
42 {
43  m_t = t;
44 }
45 
46 int32_t CDiffusionMaps::get_t() const
47 {
48  return m_t;
49 }
50 
51 const char* CDiffusionMaps::get_name() const
52 {
53  return "DiffusionMaps";
54 };
55 
57 {
58  ASSERT(features);
59  // shorthand for simplefeatures
60  SG_REF(features);
61  // compute distance matrix
63  m_kernel->init(features,features);
65  m_kernel->cleanup();
66  SG_UNREF(features);
67  return (CFeatures*)embedding;
68 }
69 
71 {
72  int32_t i,j;
73  SGMatrix<float64_t> kernel_matrix = kernel->get_kernel_matrix();
74  ASSERT(kernel_matrix.num_rows==kernel_matrix.num_cols);
75  int32_t N = kernel_matrix.num_rows;
76 
77  float64_t* p_vector = SG_CALLOC(float64_t, N);
78  for (i=0; i<N; i++)
79  {
80  for (j=0; j<N; j++)
81  {
82  p_vector[i] += kernel_matrix.matrix[j*N+i];
83  }
84  }
85 
86  float64_t* p_matrix = SG_CALLOC(float64_t, N*N);
87  cblas_dger(CblasColMajor,N,N,1.0,p_vector,1,p_vector,1,p_matrix,N);
88  for (i=0; i<N*N; i++)
89  {
90  kernel_matrix.matrix[i] /= CMath::pow(p_matrix[i], m_t);
91  }
92  SG_FREE(p_matrix);
93 
94  for (i=0; i<N; i++)
95  {
96  p_vector[i] = 0.0;
97  for (j=0; j<N; j++)
98  {
99  p_vector[i] += kernel_matrix.matrix[j*N+i];
100  }
101  p_vector[i] = CMath::sqrt(p_vector[i]);
102  }
103  float64_t ppt = cblas_ddot(N,p_vector,1,p_vector,1);
104  SG_FREE(p_vector);
105 
106  for (i=0; i<N*N; i++)
107  {
108  kernel_matrix.matrix[i] /= ppt;
109  }
110 
111 
112  float64_t* s_values = SG_MALLOC(float64_t, N);
113 
114  float64_t* kkt_matrix = SG_MALLOC(float64_t, N*N);
115 
116  cblas_dgemm(CblasColMajor,CblasTrans,CblasNoTrans,
117  N,N,N,
118  1.0,kernel_matrix.matrix,N,
119  kernel_matrix.matrix,N,
120  0.0,kkt_matrix,N);
121 
122  int32_t info = 0;
123 
124  wrap_dsyevr('V','U',N,kkt_matrix,N,N-m_target_dim,N,s_values,kernel_matrix.matrix,&info);
125  if (info)
126  SG_ERROR("DGESVD failed with %d code", info);
127 
128  SG_FREE(s_values);
129 /*
130  int32_t info = 0;
131  wrap_dgesvd('O','N',N,N,kernel_matrix.matrix,N,s_values,NULL,1,NULL,1,&info);
132 */
133 
134  SG_FREE(kkt_matrix);
135  SGMatrix<float64_t> new_feature_matrix = SGMatrix<float64_t>(m_target_dim,N);
136 
137  for (i=0; i<m_target_dim; i++)
138  {
139  for (j=0; j<N; j++)
140  new_feature_matrix[j*m_target_dim+i] = kernel_matrix.matrix[(m_target_dim-i-1)*N+j]/kernel_matrix.matrix[(m_target_dim)*N+j];
141  }
142  kernel_matrix.destroy_matrix();
143 
144  return new CSimpleFeatures<float64_t>(new_feature_matrix);
145 }
146 #endif /* HAVE_LAPACK */

SHOGUN Machine Learning Toolbox - Documentation