26 using namespace shogun;
28 #ifndef DOXYGEN_SHOULD_SKIP_THIS
29 struct LK_RECONSTRUCTION_THREAD_PARAM
42 const int32_t* neighborhood_matrix;
55 class KLLE_COVERTREE_POINT
61 this->point_index = index;
62 this->kernel_matrix = dmatrix;
63 this->kii = dmatrix[index*dmatrix.
num_rows+index];
66 inline double distance(
const KLLE_COVERTREE_POINT& p)
const
68 int32_t N = kernel_matrix.num_rows;
69 return kii+kernel_matrix[p.point_index*N+p.point_index]-2.0*kernel_matrix[point_index*N+p.point_index];
72 inline bool operator==(
const KLLE_COVERTREE_POINT& p)
const
74 return (p.point_index==this->point_index);
97 return "KernelLocallyLinearEmbedding";
112 SG_ERROR(
"Number of neighbors (%d) should be less than number of objects (%d).\n",
159 pthread_t* threads =
SG_MALLOC(pthread_t, num_threads);
160 LK_RECONSTRUCTION_THREAD_PARAM* parameters =
SG_MALLOC(LK_RECONSTRUCTION_THREAD_PARAM, num_threads);
162 pthread_attr_init(&attr);
163 pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
165 int32_t num_threads = 1;
173 for (t=0; t<num_threads; t++)
175 parameters[t].idx_start = t;
176 parameters[t].idx_step = num_threads;
177 parameters[t].idx_stop = N;
178 parameters[t].m_k =
m_k;
180 parameters[t].neighborhood_matrix = neighborhood_matrix.
matrix;
181 parameters[t].kernel_matrix = kernel_matrix.
matrix;
182 parameters[t].local_gram_matrix = local_gram_matrix+(
m_k*
m_k)*t;
183 parameters[t].id_vector = id_vector+
m_k*t;
184 parameters[t].W_matrix = W_matrix;
188 for (t=0; t<num_threads; t++)
189 pthread_join(threads[t], NULL);
190 pthread_attr_destroy(&attr);
194 LK_RECONSTRUCTION_THREAD_PARAM single_thread_param;
195 single_thread_param.idx_start = 0;
196 single_thread_param.idx_step = 1;
197 single_thread_param.idx_stop = N;
198 single_thread_param.m_k =
m_k;
199 single_thread_param.N = N;
200 single_thread_param.neighborhood_matrix = neighborhood_matrix.
matrix;
201 single_thread_param.local_gram_matrix = local_gram_matrix;
202 single_thread_param.kernel_matrix = kernel_matrix.
matrix;
203 single_thread_param.id_vector = id_vector;
204 single_thread_param.W_matrix = W_matrix;
217 LK_RECONSTRUCTION_THREAD_PARAM* parameters = (LK_RECONSTRUCTION_THREAD_PARAM*)p;
218 int32_t idx_start = parameters->idx_start;
219 int32_t idx_step = parameters->idx_step;
220 int32_t idx_stop = parameters->idx_stop;
221 int32_t
m_k = parameters->m_k;
222 int32_t N = parameters->N;
223 const int32_t* neighborhood_matrix = parameters->neighborhood_matrix;
224 float64_t* local_gram_matrix = parameters->local_gram_matrix;
225 const float64_t* kernel_matrix = parameters->kernel_matrix;
226 float64_t* id_vector = parameters->id_vector;
227 float64_t* W_matrix = parameters->W_matrix;
228 float64_t reconstruction_shift = parameters->reconstruction_shift;
233 for (i=idx_start; i<idx_stop; i+=idx_step)
235 for (j=0; j<
m_k; j++)
237 for (k=0; k<
m_k; k++)
238 local_gram_matrix[j*m_k+k] =
239 kernel_matrix[i*N+i] -
240 kernel_matrix[i*N+neighborhood_matrix[i*m_k+j]] -
241 kernel_matrix[i*N+neighborhood_matrix[i*m_k+k]] +
242 kernel_matrix[neighborhood_matrix[i*m_k+j]*N+neighborhood_matrix[i*m_k+k]];
245 for (j=0; j<
m_k; j++)
250 for (j=0; j<
m_k; j++)
251 trace += local_gram_matrix[j*m_k+j];
254 for (j=0; j<
m_k; j++)
255 local_gram_matrix[j*m_k+j] += reconstruction_shift*trace;
257 clapack_dposv(CblasColMajor,CblasLower,m_k,1,local_gram_matrix,m_k,id_vector,m_k);
261 for (j=0; j<
m_k; j++)
262 norming += id_vector[j];
264 cblas_dscal(m_k,1.0/norming,id_vector,1);
266 memset(local_gram_matrix,0,
sizeof(
float64_t)*m_k*m_k);
267 cblas_dger(CblasColMajor,m_k,m_k,1.0,id_vector,1,id_vector,1,local_gram_matrix,m_k);
270 W_matrix[N*i+i] += 1.0;
271 for (j=0; j<
m_k; j++)
273 W_matrix[N*i+neighborhood_matrix[i*m_k+j]] -= id_vector[j];
274 W_matrix[N*neighborhood_matrix[i*m_k+j]+i] -= id_vector[j];
276 for (j=0; j<
m_k; j++)
278 for (k=0; k<
m_k; k++)
279 W_matrix[N*neighborhood_matrix[i*m_k+j]+neighborhood_matrix[i*m_k+k]]+=local_gram_matrix[j*m_k+k];
290 int32_t* neighborhood_matrix =
SG_MALLOC(int32_t, N*k);
294 max_dist =
CMath::max(max_dist,kernel_matrix[i*N+i]);
296 std::vector<KLLE_COVERTREE_POINT> vectors;
299 vectors.push_back(KLLE_COVERTREE_POINT(i,kernel_matrix));
305 std::vector<KLLE_COVERTREE_POINT> neighbors =
308 ASSERT(neighbors.size()>=unsigned(k+1));
310 for (std::size_t m=1; m<unsigned(k+1); m++)
311 neighborhood_matrix[i*k+m-1] = neighbors[m].point_index;