00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #include <shogun/converter/LocalityPreservingProjections.h>
00012 #ifdef HAVE_LAPACK
00013 #include <shogun/mathematics/arpack.h>
00014 #include <shogun/mathematics/lapack.h>
00015 #include <shogun/lib/FibonacciHeap.h>
00016 #include <shogun/mathematics/Math.h>
00017 #include <shogun/io/SGIO.h>
00018 #include <shogun/distance/Distance.h>
00019 #include <shogun/lib/Signal.h>
00020
00021 using namespace shogun;
00022
00023 CLocalityPreservingProjections::CLocalityPreservingProjections() :
00024 CLaplacianEigenmaps()
00025 {
00026 }
00027
00028 CLocalityPreservingProjections::~CLocalityPreservingProjections()
00029 {
00030 }
00031
00032 const char* CLocalityPreservingProjections::get_name() const
00033 {
00034 return "LocalityPreservingProjections";
00035 };
00036
00037 CDenseFeatures<float64_t>* CLocalityPreservingProjections::construct_embedding(CFeatures* features,
00038 SGMatrix<float64_t> W_matrix)
00039 {
00040 CDenseFeatures<float64_t>* simple_features = (CDenseFeatures<float64_t>*)features;
00041 ASSERT(simple_features);
00042 int i,j;
00043 int N = simple_features->get_num_vectors();
00044 int dim = simple_features->get_num_features();
00045
00046 float64_t* D_diag_vector = SG_CALLOC(float64_t, N);
00047 for (i=0; i<N; i++)
00048 {
00049 for (j=0; j<N; j++)
00050 D_diag_vector[i] += W_matrix[i*N+j];
00051 }
00052
00053
00054 for (i=0; i<N*N; i++)
00055 if (W_matrix[i]>0.0)
00056 W_matrix[i] *= -1.0;
00057
00058 for (i=0; i<N; i++)
00059 W_matrix[i*N+i] += D_diag_vector[i];
00060
00061 SGMatrix<float64_t> feature_matrix = simple_features->get_feature_matrix();
00062 float64_t* XTM = SG_MALLOC(float64_t, dim*N);
00063 float64_t* lhs_M = SG_MALLOC(float64_t, dim*dim);
00064 float64_t* rhs_M = SG_MALLOC(float64_t, dim*dim);
00065
00066 cblas_dgemm(CblasColMajor,CblasNoTrans,CblasNoTrans,dim,N,N,1.0,feature_matrix.matrix,dim,W_matrix.matrix,N,0.0,XTM,dim);
00067 cblas_dgemm(CblasColMajor,CblasNoTrans,CblasTrans,dim,dim,N,1.0,XTM,dim,feature_matrix.matrix,dim,0.0,lhs_M,dim);
00068
00069 for (i=0; i<N; i++)
00070 cblas_dscal(dim,CMath::sqrt(D_diag_vector[i]),feature_matrix.matrix+i*dim,1);
00071
00072 cblas_dgemm(CblasColMajor,CblasNoTrans,CblasTrans,dim,dim,N,1.0,feature_matrix.matrix,dim,feature_matrix.matrix,dim,0.0,rhs_M,dim);
00073
00074 for (i=0; i<N; i++)
00075 cblas_dscal(dim,1.0/D_diag_vector[i],feature_matrix.matrix+i*dim,1);
00076
00077 float64_t* evals = SG_MALLOC(float64_t, dim);
00078 float64_t* evectors = SG_MALLOC(float64_t, m_target_dim*dim);
00079 int32_t info = 0;
00080 #ifdef HAVE_ARPACK
00081 arpack_dsxupd(lhs_M,rhs_M,false,dim,m_target_dim,"LA",false,3,true,false,-1e-9,0.0,
00082 evals,evectors,info);
00083 #else
00084 wrap_dsygvx(1,'V','U',dim,lhs_M,dim,rhs_M,dim,1,m_target_dim,evals,evectors,&info);
00085 #endif
00086 SG_FREE(lhs_M);
00087 SG_FREE(rhs_M);
00088
00089 SG_FREE(evals);
00090
00091 if (info!=0)
00092 SG_ERROR("Failed to solve eigenproblem (%d)\n",info);
00093
00094 cblas_dgemm(CblasColMajor,CblasTrans,CblasNoTrans,N,m_target_dim,dim,1.0,feature_matrix.matrix,dim,evectors,dim,0.0,XTM,N);
00095 SG_FREE(evectors);
00096
00097 SGMatrix<float64_t> new_features(m_target_dim,N);
00098 for (i=0; i<m_target_dim; i++)
00099 {
00100 for (j=0; j<N; j++)
00101 new_features[j*m_target_dim+i] = XTM[i*N+j];
00102 }
00103 SG_FREE(D_diag_vector);
00104 SG_FREE(XTM);
00105 return new CDenseFeatures<float64_t>(new_features);
00106 }
00107
00108 #endif