SHOGUN
v2.0.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 Sergey Lisitsyn 00009 */ 00010 00011 #include <shogun/converter/LinearLocalTangentSpaceAlignment.h> 00012 #ifdef HAVE_LAPACK 00013 #include <shogun/mathematics/arpack.h> 00014 #include <shogun/mathematics/lapack.h> 00015 #include <shogun/mathematics/Math.h> 00016 #include <shogun/io/SGIO.h> 00017 #include <shogun/lib/Time.h> 00018 #include <shogun/distance/Distance.h> 00019 #include <shogun/lib/Signal.h> 00020 00021 using namespace shogun; 00022 00023 CLinearLocalTangentSpaceAlignment::CLinearLocalTangentSpaceAlignment() : 00024 CLocalTangentSpaceAlignment() 00025 { 00026 } 00027 00028 CLinearLocalTangentSpaceAlignment::~CLinearLocalTangentSpaceAlignment() 00029 { 00030 } 00031 00032 const char* CLinearLocalTangentSpaceAlignment::get_name() const 00033 { 00034 return "LinearLocalTangentSpaceAlignment"; 00035 } 00036 00037 SGMatrix<float64_t> CLinearLocalTangentSpaceAlignment::construct_embedding(CFeatures* features, SGMatrix<float64_t> matrix, int dimension) 00038 { 00039 CDenseFeatures<float64_t>* simple_features = (CDenseFeatures<float64_t>*)features; 00040 ASSERT(simple_features); 00041 int i,j; 00042 00043 SGMatrix<float64_t> feature_matrix = simple_features->get_feature_matrix().clone(); 00044 int N=feature_matrix.num_cols; 00045 int dim=feature_matrix.num_rows; 00046 ASSERT(dimension<=dim); 00047 float64_t* XTM = SG_MALLOC(float64_t, dim*N); 00048 float64_t* lhs_M = SG_MALLOC(float64_t, dim*dim); 00049 float64_t* rhs_M = SG_MALLOC(float64_t, dim*dim); 00050 SGMatrix<float64_t>::center_matrix(matrix.matrix,N,N); 00051 00052 cblas_dgemm(CblasColMajor,CblasNoTrans,CblasNoTrans,dim,N,N,1.0,feature_matrix.matrix,dim,matrix.matrix,N,0.0,XTM,dim); 00053 cblas_dgemm(CblasColMajor,CblasNoTrans,CblasTrans,dim,dim,N,1.0,XTM,dim,feature_matrix.matrix,dim,0.0,lhs_M,dim); 00054 00055 float64_t* mean_vector = SG_CALLOC(float64_t, dim); 00056 for (i=0; i<N; i++) 00057 cblas_daxpy(dim,1.0,feature_matrix.matrix+i*dim,1,mean_vector,1); 00058 00059 cblas_dscal(dim,1.0/N,mean_vector,1); 00060 00061 for (i=0; i<N; i++) 00062 cblas_daxpy(dim,-1.0,mean_vector,1,feature_matrix.matrix+i*dim,1); 00063 00064 SG_FREE(mean_vector); 00065 00066 cblas_dgemm(CblasColMajor,CblasNoTrans,CblasTrans,dim,dim,N,1.0,feature_matrix.matrix,dim,feature_matrix.matrix,dim,0.0,rhs_M,dim); 00067 00068 for (i=0; i<dim; i++) rhs_M[i*dim+i] += 1e-6; 00069 00070 float64_t* evals = SG_MALLOC(float64_t, dim); 00071 float64_t* evectors = SG_MALLOC(float64_t, dimension*dim); 00072 int32_t info = 0; 00073 #ifdef HAVE_ARPACK 00074 arpack_dsxupd(lhs_M,rhs_M,false,dim,dimension,"LA",false,3,true,false,m_nullspace_shift,0.0, 00075 evals,evectors,info); 00076 #else 00077 wrap_dsygvx(1,'V','U',dim,lhs_M,dim,rhs_M,dim,dim-dimension+1,dim,evals,evectors,&info); 00078 #endif 00079 SG_FREE(lhs_M); 00080 SG_FREE(rhs_M); 00081 SG_FREE(evals); 00082 if (info!=0) SG_ERROR("Failed to solve eigenproblem (%d)\n",info); 00083 00084 for (i=0; i<dimension/2; i++) 00085 { 00086 cblas_dswap(dim,evectors+i*dim,1,evectors+(dimension-i-1)*dim,1); 00087 } 00088 00089 cblas_dgemm(CblasColMajor,CblasTrans,CblasNoTrans,N,dimension,dim,1.0,feature_matrix.matrix,dim,evectors,dim,0.0,XTM,N); 00090 SG_FREE(evectors); 00091 00092 float64_t* new_features = SG_MALLOC(float64_t, dimension*N); 00093 for (i=0; i<dimension; i++) 00094 { 00095 for (j=0; j<N; j++) 00096 new_features[j*dimension+i] = XTM[i*N+j]; 00097 } 00098 SG_FREE(XTM); 00099 return SGMatrix<float64_t>(new_features,dimension,N); 00100 } 00101 00102 #endif /* HAVE_LAPACK */