SHOGUN  v2.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups 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
19 #include <shogun/io/SGIO.h>
20 #include <shogun/kernel/Kernel.h>
21 #include <shogun/lib/Signal.h>
22 #include <shogun/lib/Time.h>
23 
24 using namespace shogun;
25 
28 {
29  m_t = 10;
30  set_kernel(new CGaussianKernel(10,1.0));
31 
32  init();
33 }
34 
36 {
37  SG_ADD(&m_t, "t", "number of steps", MS_AVAILABLE);
38 }
39 
41 {
42 }
43 
44 void CDiffusionMaps::set_t(int32_t t)
45 {
46  m_t = t;
47 }
48 
49 int32_t CDiffusionMaps::get_t() const
50 {
51  return m_t;
52 }
53 
54 const char* CDiffusionMaps::get_name() const
55 {
56  return "DiffusionMaps";
57 };
58 
60 {
61  ASSERT(features);
62  // shorthand for simplefeatures
63  SG_REF(features);
64  // compute distance matrix
66  m_kernel->init(features,features);
68  m_kernel->cleanup();
69  SG_UNREF(features);
70  return (CFeatures*)embedding;
71 }
72 
74 {
75 #ifdef HAVE_ARPACK
76  bool use_arpack = true;
77 #else
78  bool use_arpack = false;
79 #endif
80  int32_t i,j;
81  SGMatrix<float64_t> kernel_matrix = kernel->get_kernel_matrix();
82  ASSERT(kernel_matrix.num_rows==kernel_matrix.num_cols);
83  int32_t N = kernel_matrix.num_rows;
84 
85  float64_t* p_vector = SG_CALLOC(float64_t, N);
86  for (i=0; i<N; i++)
87  {
88  for (j=0; j<N; j++)
89  {
90  p_vector[i] += kernel_matrix.matrix[j*N+i];
91  }
92  }
93  //CMath::display_matrix(kernel_matrix.matrix,N,N,"K");
94  for (i=0; i<N; i++)
95  {
96  for (j=0; j<N; j++)
97  {
98  kernel_matrix.matrix[i*N+j] /= CMath::pow(p_vector[i]*p_vector[j], m_t);
99  }
100  }
101  //CMath::display_matrix(kernel_matrix.matrix,N,N,"K");
102 
103  for (i=0; i<N; i++)
104  {
105  p_vector[i] = 0.0;
106  for (j=0; j<N; j++)
107  {
108  p_vector[i] += kernel_matrix.matrix[j*N+i];
109  }
110  p_vector[i] = CMath::sqrt(p_vector[i]);
111  }
112 
113  for (i=0; i<N; i++)
114  {
115  for (j=0; j<N; j++)
116  {
117  kernel_matrix.matrix[i*N+j] /= p_vector[i]*p_vector[j];
118  }
119  }
120 
121  SG_FREE(p_vector);
122  float64_t* s_values = SG_MALLOC(float64_t, N);
123 
124  int32_t info = 0;
125  SGMatrix<float64_t> new_feature_matrix = SGMatrix<float64_t>(m_target_dim,N);
126 
127  if (use_arpack)
128  {
129 #ifdef HAVE_ARPACK
130  arpack_dsxupd(kernel_matrix.matrix,NULL,false,N,m_target_dim,"LA",false,1,false,true,0.0,0.0,s_values,kernel_matrix.matrix,info);
131 #endif /* HAVE_ARPACK */
132  for (i=0; i<m_target_dim; i++)
133  {
134  for (j=0; j<N; j++)
135  new_feature_matrix[j*m_target_dim+i] = kernel_matrix[j*m_target_dim+i];
136  }
137  }
138  else
139  {
140  SG_WARNING("LAPACK does not provide efficient routines to construct embedding (this may take time). Consider installing ARPACK.");
141  wrap_dgesvd('O','N',N,N,kernel_matrix.matrix,N,s_values,NULL,1,NULL,1,&info);
142  for (i=0; i<m_target_dim; i++)
143  {
144  for (j=0; j<N; j++)
145  new_feature_matrix[j*m_target_dim+i] =
146  kernel_matrix[(m_target_dim-i-1)*N+j];
147  }
148  }
149  if (info)
150  SG_ERROR("Eigenproblem solving failed with %d code", info);
151 
152  SG_FREE(s_values);
153 
154  return new CDenseFeatures<float64_t>(new_feature_matrix);
155 }
156 #endif /* HAVE_LAPACK */

SHOGUN Machine Learning Toolbox - Documentation