SHOGUN
v1.1.0
|
00001 /* 00002 * This program is free software; you can redistribute it and/or modify 00003 * it under the terms of the GNU General Public License as published by 00004 * the Free Software Foundation; either version 3 of the License, or 00005 * (at your option) any later version. 00006 * 00007 * Written (W) 2011 Sergey Lisitsyn 00008 * Copyright (C) 2011 Berlin Institute of Technology and Max-Planck-Society 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 CSimpleFeatures<float64_t>* CLocalityPreservingProjections::construct_embedding(CFeatures* features, 00038 SGMatrix<float64_t> W_matrix) 00039 { 00040 CSimpleFeatures<float64_t>* simple_features = (CSimpleFeatures<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 // W = -W 00054 for (i=0; i<N*N; i++) 00055 if (W_matrix[i]>0.0) 00056 W_matrix[i] *= -1.0; 00057 // W = W + D 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,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,-1e-9,0.0, 00082 evals,evectors,info); 00083 #else 00084 wrap_dsygvx(1,'V','U',dim,lhs_M,dim,rhs_M,dim,dim-m_target_dim+1,dim,evals,evectors,&info); 00085 #endif 00086 SG_FREE(lhs_M); 00087 SG_FREE(rhs_M); 00088 SG_FREE(evals); 00089 if (info!=0) SG_ERROR("Failed to solve eigenproblem (%d)\n",info); 00090 00091 cblas_dgemm(CblasColMajor,CblasTrans,CblasNoTrans,N,m_target_dim,dim,1.0,feature_matrix.matrix,dim,evectors,dim,0.0,XTM,N); 00092 SG_FREE(evectors); 00093 00094 SGMatrix<float64_t> new_features(m_target_dim,N); 00095 for (i=0; i<m_target_dim; i++) 00096 { 00097 for (j=0; j<N; j++) 00098 new_features[j*m_target_dim+i] = XTM[i*N+j]; 00099 } 00100 SG_FREE(D_diag_vector); 00101 SG_FREE(XTM); 00102 return new CSimpleFeatures<float64_t>(new_features); 00103 } 00104 00105 #endif /* HAVE_LAPACK */