00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #include <shogun/converter/DiffusionMaps.h>
00012 #include <shogun/converter/EmbeddingConverter.h>
00013 #include <shogun/lib/config.h>
00014 #ifdef HAVE_LAPACK
00015 #include <shogun/mathematics/lapack.h>
00016 #include <shogun/mathematics/Math.h>
00017 #include <shogun/io/SGIO.h>
00018 #include <shogun/kernel/Kernel.h>
00019 #include <shogun/lib/Signal.h>
00020 #include <shogun/lib/Time.h>
00021
00022 using namespace shogun;
00023
00024 CDiffusionMaps::CDiffusionMaps() :
00025 CEmbeddingConverter()
00026 {
00027 m_t = 10;
00028
00029 init();
00030 }
00031
00032 void CDiffusionMaps::init()
00033 {
00034 m_parameters->add(&m_t, "t", "number of steps");
00035 }
00036
00037 CDiffusionMaps::~CDiffusionMaps()
00038 {
00039 }
00040
00041 void CDiffusionMaps::set_t(int32_t t)
00042 {
00043 m_t = t;
00044 }
00045
00046 int32_t CDiffusionMaps::get_t() const
00047 {
00048 return m_t;
00049 }
00050
00051 const char* CDiffusionMaps::get_name() const
00052 {
00053 return "DiffusionMaps";
00054 };
00055
00056 CFeatures* CDiffusionMaps::apply(CFeatures* features)
00057 {
00058 ASSERT(features);
00059
00060 SG_REF(features);
00061
00062 ASSERT(m_kernel);
00063 m_kernel->init(features,features);
00064 CSimpleFeatures<float64_t>* embedding = embed_kernel(m_kernel);
00065 m_kernel->cleanup();
00066 SG_UNREF(features);
00067 return (CFeatures*)embedding;
00068 }
00069
00070 CSimpleFeatures<float64_t>* CDiffusionMaps::embed_kernel(CKernel* kernel)
00071 {
00072 int32_t i,j;
00073 SGMatrix<float64_t> kernel_matrix = kernel->get_kernel_matrix();
00074 ASSERT(kernel_matrix.num_rows==kernel_matrix.num_cols);
00075 int32_t N = kernel_matrix.num_rows;
00076
00077 float64_t* p_vector = SG_CALLOC(float64_t, N);
00078 for (i=0; i<N; i++)
00079 {
00080 for (j=0; j<N; j++)
00081 {
00082 p_vector[i] += kernel_matrix.matrix[j*N+i];
00083 }
00084 }
00085
00086 float64_t* p_matrix = SG_CALLOC(float64_t, N*N);
00087 cblas_dger(CblasColMajor,N,N,1.0,p_vector,1,p_vector,1,p_matrix,N);
00088 for (i=0; i<N*N; i++)
00089 {
00090 kernel_matrix.matrix[i] /= CMath::pow(p_matrix[i], m_t);
00091 }
00092 SG_FREE(p_matrix);
00093
00094 for (i=0; i<N; i++)
00095 {
00096 p_vector[i] = 0.0;
00097 for (j=0; j<N; j++)
00098 {
00099 p_vector[i] += kernel_matrix.matrix[j*N+i];
00100 }
00101 p_vector[i] = CMath::sqrt(p_vector[i]);
00102 }
00103 float64_t ppt = cblas_ddot(N,p_vector,1,p_vector,1);
00104 SG_FREE(p_vector);
00105
00106 for (i=0; i<N*N; i++)
00107 {
00108 kernel_matrix.matrix[i] /= ppt;
00109 }
00110
00111
00112 float64_t* s_values = SG_MALLOC(float64_t, N);
00113
00114 float64_t* kkt_matrix = SG_MALLOC(float64_t, N*N);
00115
00116 cblas_dgemm(CblasColMajor,CblasTrans,CblasNoTrans,
00117 N,N,N,
00118 1.0,kernel_matrix.matrix,N,
00119 kernel_matrix.matrix,N,
00120 0.0,kkt_matrix,N);
00121
00122 int32_t info = 0;
00123
00124 wrap_dsyevr('V','U',N,kkt_matrix,N,N-m_target_dim,N,s_values,kernel_matrix.matrix,&info);
00125 if (info)
00126 SG_ERROR("DGESVD failed with %d code", info);
00127
00128 SG_FREE(s_values);
00129
00130
00131
00132
00133
00134 SG_FREE(kkt_matrix);
00135 SGMatrix<float64_t> new_feature_matrix = SGMatrix<float64_t>(m_target_dim,N);
00136
00137 for (i=0; i<m_target_dim; i++)
00138 {
00139 for (j=0; j<N; j++)
00140 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];
00141 }
00142 kernel_matrix.destroy_matrix();
00143
00144 return new CSimpleFeatures<float64_t>(new_feature_matrix);
00145 }
00146 #endif