Skip to content

Commit

Permalink
Browse files Browse the repository at this point in the history
CoverTree for EDRT
  • Loading branch information
lisitsyn committed Jan 30, 2012
1 parent 689b5ed commit a253451
Show file tree
Hide file tree
Showing 6 changed files with 760 additions and 230 deletions.
59 changes: 41 additions & 18 deletions src/shogun/converter/Isomap.cpp
Expand Up @@ -16,6 +16,7 @@
#include <shogun/io/SGIO.h>
#include <shogun/base/Parallel.h>
#include <shogun/lib/Signal.h>
#include <shogun/lib/CoverTree.h>

#ifdef HAVE_PTHREAD
#include <pthread.h>
Expand Down Expand Up @@ -51,6 +52,31 @@ struct DIJKSTRA_THREAD_PARAM
/// (f)rontier bool array
bool* f;
};

class ISOMAP_COVERTREE_POINT
{
public:

ISOMAP_COVERTREE_POINT(int32_t index, const SGMatrix<float64_t>& dmatrix)
{
point_index = index;
distance_matrix = dmatrix;
}

inline double distance(const ISOMAP_COVERTREE_POINT& p) const
{
return distance_matrix[point_index*distance_matrix.num_rows+p.point_index];
}

inline bool operator==(const ISOMAP_COVERTREE_POINT& p) const
{
return (p.point_index==point_index);
}

int32_t point_index;
SGMatrix<float64_t> distance_matrix;
};

#endif /* DOXYGEN_SHOULD_SKIP_THIS */

CIsomap::CIsomap() : CMultidimensionalScaling()
Expand Down Expand Up @@ -109,29 +135,26 @@ SGMatrix<float64_t> CIsomap::isomap_distance(SGMatrix<float64_t> D_matrix)
// cut by k-nearest neighbors
int32_t* edges_idx_matrix = SG_MALLOC(int32_t, N*m_k);
float64_t* edges_matrix = SG_MALLOC(float64_t, N*m_k);

// query neighbors and edges to neighbors
CFibonacciHeap* heap = new CFibonacciHeap(N);
for (i=0; i<N; i++)
{
// insert distances to heap
for (j=0; j<N; j++)
heap->insert(j,D_matrix[i*N+j]);

// extract nearest neighbor: the jth object itself
heap->extract_min(tmp);
float64_t max_dist = CMath::max(D_matrix.matrix,N*N);
CoverTree<ISOMAP_COVERTREE_POINT>* coverTree = new CoverTree<ISOMAP_COVERTREE_POINT>(max_dist);

// extract m_k neighbors and distances
for (j=0; j<m_k; j++)
for (i=0; i<N; i++)
coverTree->insert(ISOMAP_COVERTREE_POINT(i,D_matrix));

for (i=0; i<N; i++)
{
ISOMAP_COVERTREE_POINT origin(i,D_matrix);
std::vector<ISOMAP_COVERTREE_POINT> neighbors =
coverTree->kNearestNeighbors(origin,m_k+1);
for (std::size_t m=1; m<neighbors.size(); m++)
{
edges_idx_matrix[i*m_k+j] = heap->extract_min(tmp);
edges_matrix[i*m_k+j] = tmp;
edges_matrix[i*m_k+m-1] = origin.distance(neighbors[m]);
edges_idx_matrix[i*m_k+m-1] = neighbors[m].point_index;
}
// clear heap
heap->clear();
}
// cleanup
delete heap;

delete coverTree;

#ifdef HAVE_PTHREAD

Expand Down
145 changes: 39 additions & 106 deletions src/shogun/converter/KernelLocallyLinearEmbedding.cpp
Expand Up @@ -52,27 +52,33 @@ struct LK_RECONSTRUCTION_THREAD_PARAM
float64_t* W_matrix;
};

struct K_NEIGHBORHOOD_THREAD_PARAM
class KLLE_COVERTREE_POINT
{
/// starting index of loop
int32_t idx_start;
/// step of loop
int32_t idx_step;
/// end index of loop
int32_t idx_stop;
/// number of vectors
int32_t N;
/// number of neighbors
int32_t m_k;
/// fibonacci heaps
CFibonacciHeap* heap;
/// kernel matrix
const float64_t* kernel_matrix;
/// kernel diag
const float64_t* kernel_diag;
/// matrix containing neighbors indexes
int32_t* neighborhood_matrix;
public:

KLLE_COVERTREE_POINT(int32_t index, const SGMatrix<float64_t>& dmatrix)
{
this->point_index = index;
this->kernel_matrix = dmatrix;
this->kii = dmatrix[index*dmatrix.num_rows+index];
}

inline double distance(const KLLE_COVERTREE_POINT& p) const
{
return kernel_matrix[p.point_index*kernel_matrix.num_rows+p.point_index]+kii-
2.0*kernel_matrix[point_index*kernel_matrix.num_rows+p.point_index];
}

inline bool operator==(const KLLE_COVERTREE_POINT& p) const
{
return (p.point_index==this->point_index);
}

int32_t point_index;
float64_t kii;
SGMatrix<float64_t> kernel_matrix;
};

#endif /* DOXYGEN_SHOULD_SKIP_THIS */

CKernelLocallyLinearEmbedding::CKernelLocallyLinearEmbedding() :
Expand Down Expand Up @@ -281,102 +287,29 @@ void* CKernelLocallyLinearEmbedding::run_linearreconstruction_thread(void* p)

SGMatrix<int32_t> CKernelLocallyLinearEmbedding::get_neighborhood_matrix(SGMatrix<float64_t> kernel_matrix, int32_t k)
{
int32_t t;
int32_t i;
int32_t N = kernel_matrix.num_cols;
// init matrix and heap to be used

int32_t* neighborhood_matrix = SG_MALLOC(int32_t, N*k);
float64_t* kernel_diag = SG_MALLOC(float64_t, N);
// dirty
for (t=0; t<N; t++)
kernel_diag[t] = kernel_matrix.matrix[t*N+t];

float64_t max_dist = CMath::max(kernel_matrix.matrix,N*N);

#ifdef HAVE_PTHREAD
int32_t num_threads = parallel->get_num_threads();
ASSERT(num_threads>0);
K_NEIGHBORHOOD_THREAD_PARAM* parameters = SG_MALLOC(K_NEIGHBORHOOD_THREAD_PARAM, num_threads);
pthread_t* threads = SG_MALLOC(pthread_t, num_threads);
pthread_attr_t attr;
pthread_attr_init(&attr);
pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
#else
int32_t num_threads = 1;
#endif
CFibonacciHeap** heaps = SG_MALLOC(CFibonacciHeap*, num_threads);
for (t=0; t<num_threads; t++)
heaps[t] = new CFibonacciHeap(N);
CoverTree<KLLE_COVERTREE_POINT>* coverTree = new CoverTree<KLLE_COVERTREE_POINT>(max_dist);

#ifdef HAVE_PTHREAD
for (t=0; t<num_threads; t++)
for (i=0; i<N; i++)
coverTree->insert(KLLE_COVERTREE_POINT(i,kernel_matrix));

for (i=0; i<N; i++)
{
parameters[t].idx_start = t;
parameters[t].idx_step = num_threads;
parameters[t].idx_stop = N;
parameters[t].m_k = k;
parameters[t].N = N;
parameters[t].heap = heaps[t];
parameters[t].neighborhood_matrix = neighborhood_matrix;
parameters[t].kernel_matrix = kernel_matrix.matrix;
parameters[t].kernel_diag = kernel_diag;
pthread_create(&threads[t], &attr, run_neighborhood_thread, (void*)&parameters[t]);
std::vector<KLLE_COVERTREE_POINT> neighbors =
coverTree->kNearestNeighbors(KLLE_COVERTREE_POINT(i,kernel_matrix),k+1);
for (std::size_t m=1; m<neighbors.size(); m++)
neighborhood_matrix[i*k+m-1] = neighbors[m].point_index;
}
for (t=0; t<num_threads; t++)
pthread_join(threads[t], NULL);
pthread_attr_destroy(&attr);
SG_FREE(threads);
SG_FREE(parameters);
#else
K_NEIGHBORHOOD_THREAD_PARAM single_thread_param;
single_thread_param.idx_start = 0;
single_thread_param.idx_step = 1;
single_thread_param.idx_stop = N;
single_thread_param.m_k = k;
single_thread_param.N = N;
single_thread_param.heap = heaps[0]
single_thread_param.neighborhood_matrix = neighborhood_matrix;
single_thread_param.kernel_matrix = kernel_matrix.matrix;
single_thread_param.kernel_diag = kernel_diag;
run_neighborhood_thread((void*)&single_thread_param);
#endif

SG_FREE(kernel_diag);
for (t=0; t<num_threads; t++)
delete heaps[t];
SG_FREE(heaps);
delete coverTree;

return SGMatrix<int32_t>(neighborhood_matrix,k,N);
}

void* CKernelLocallyLinearEmbedding::run_neighborhood_thread(void* p)
{
K_NEIGHBORHOOD_THREAD_PARAM* parameters = (K_NEIGHBORHOOD_THREAD_PARAM*)p;
int32_t idx_start = parameters->idx_start;
int32_t idx_step = parameters->idx_step;
int32_t idx_stop = parameters->idx_stop;
int32_t N = parameters->N;
int32_t m_k = parameters->m_k;
CFibonacciHeap* heap = parameters->heap;
const float64_t* kernel_matrix = parameters->kernel_matrix;
const float64_t* kernel_diag = parameters->kernel_diag;
int32_t* neighborhood_matrix = parameters->neighborhood_matrix;

int32_t i,j;
float64_t tmp,k_diag_i,k_diag_j;
for (i=idx_start; i<idx_stop; i+=idx_step)
{
k_diag_i = kernel_diag[i];
for (j=0; j<N; j++)
{
heap->insert(j,k_diag_i+kernel_diag[j]-2.0*kernel_matrix[i*N+j]);
}

heap->extract_min(tmp);

for (j=0; j<m_k; j++)
neighborhood_matrix[i*m_k+j] = heap->extract_min(tmp);

heap->clear();
}

return NULL;
}
#endif /* HAVE_LAPACK */
6 changes: 1 addition & 5 deletions src/shogun/converter/KernelLocallyLinearEmbedding.h
Expand Up @@ -16,6 +16,7 @@
#include <shogun/features/Features.h>
#include <shogun/distance/Distance.h>
#include <shogun/kernel/Kernel.h>
#include <shogun/lib/CoverTree.h>

namespace shogun
{
Expand Down Expand Up @@ -85,11 +86,6 @@ class CKernelLocallyLinearEmbedding: public CLocallyLinearEmbedding
/// THREADS
protected:

/** runs neighborhood determination thread
* @param p thread params
*/
static void* run_neighborhood_thread(void* p);

/** runs linear reconstruction thread
* @param p thread params
*/
Expand Down

0 comments on commit a253451

Please sign in to comment.