25 using namespace shogun;
27 #ifndef DOXYGEN_SHOULD_SKIP_THIS
28 struct LK_RECONSTRUCTION_THREAD_PARAM
41 const int32_t* neighborhood_matrix;
54 struct K_NEIGHBORHOOD_THREAD_PARAM
71 int32_t* neighborhood_matrix;
88 return "KernelLocallyLinearEmbedding";
103 SG_ERROR(
"Number of neighbors (%d) should be less than number of objects (%d).\n",
139 pthread_t* threads =
SG_MALLOC(pthread_t, num_threads);
140 LK_RECONSTRUCTION_THREAD_PARAM* parameters =
SG_MALLOC(LK_RECONSTRUCTION_THREAD_PARAM, num_threads);
142 pthread_attr_init(&attr);
143 pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
145 int32_t num_threads = 1;
153 for (t=0; t<num_threads; t++)
155 parameters[t].idx_start = t;
156 parameters[t].idx_step = num_threads;
157 parameters[t].idx_stop = N;
158 parameters[t].m_k =
m_k;
160 parameters[t].neighborhood_matrix = neighborhood_matrix.
matrix;
161 parameters[t].kernel_matrix = kernel_matrix.
matrix;
162 parameters[t].local_gram_matrix = local_gram_matrix+(
m_k*
m_k)*t;
163 parameters[t].id_vector = id_vector+
m_k*t;
164 parameters[t].W_matrix = W_matrix;
168 for (t=0; t<num_threads; t++)
169 pthread_join(threads[t], NULL);
170 pthread_attr_destroy(&attr);
174 LK_RECONSTRUCTION_THREAD_PARAM single_thread_param;
175 single_thread_param.idx_start = 0;
176 single_thread_param.idx_step = 1;
177 single_thread_param.idx_stop = N;
178 single_thread_param.m_k =
m_k;
179 single_thread_param.N = N;
180 single_thread_param.neighborhood_matrix = neighborhood_matrix.
matrix;
181 single_thread_param.local_gram_matrix = local_gram_matrix;
182 single_thread_param.kernel_matrix = kernel_matrix.
matrix;
183 single_thread_param.id_vector = id_vector;
184 single_thread_param.W_matrix = W_matrix;
197 LK_RECONSTRUCTION_THREAD_PARAM* parameters = (LK_RECONSTRUCTION_THREAD_PARAM*)p;
198 int32_t idx_start = parameters->idx_start;
199 int32_t idx_step = parameters->idx_step;
200 int32_t idx_stop = parameters->idx_stop;
201 int32_t
m_k = parameters->m_k;
202 int32_t N = parameters->N;
203 const int32_t* neighborhood_matrix = parameters->neighborhood_matrix;
204 float64_t* local_gram_matrix = parameters->local_gram_matrix;
205 const float64_t* kernel_matrix = parameters->kernel_matrix;
206 float64_t* id_vector = parameters->id_vector;
207 float64_t* W_matrix = parameters->W_matrix;
208 float64_t reconstruction_shift = parameters->reconstruction_shift;
213 for (i=idx_start; i<idx_stop; i+=idx_step)
215 for (j=0; j<
m_k; j++)
217 for (k=0; k<
m_k; k++)
218 local_gram_matrix[j*m_k+k] =
219 kernel_matrix[i*N+i] -
220 kernel_matrix[i*N+neighborhood_matrix[j*N+i]] -
221 kernel_matrix[i*N+neighborhood_matrix[k*N+i]] +
222 kernel_matrix[neighborhood_matrix[j*N+i]*N+neighborhood_matrix[k*N+i]];
225 for (j=0; j<
m_k; j++)
230 for (j=0; j<
m_k; j++)
231 trace += local_gram_matrix[j*m_k+j];
234 for (j=0; j<
m_k; j++)
235 local_gram_matrix[j*m_k+j] += reconstruction_shift*trace;
237 clapack_dposv(CblasColMajor,CblasLower,m_k,1,local_gram_matrix,m_k,id_vector,m_k);
241 for (j=0; j<
m_k; j++)
242 norming += id_vector[j];
244 cblas_dscal(m_k,1.0/norming,id_vector,1);
246 memset(local_gram_matrix,0,
sizeof(
float64_t)*m_k*m_k);
247 cblas_dger(CblasColMajor,m_k,m_k,1.0,id_vector,1,id_vector,1,local_gram_matrix,m_k);
250 W_matrix[N*i+i] += 1.0;
251 for (j=0; j<
m_k; j++)
253 W_matrix[N*i+neighborhood_matrix[j*N+i]] -= id_vector[j];
254 W_matrix[N*neighborhood_matrix[j*N+i]+i] -= id_vector[j];
256 for (j=0; j<
m_k; j++)
258 for (k=0; k<
m_k; k++)
259 W_matrix[N*neighborhood_matrix[j*N+i]+neighborhood_matrix[k*N+i]]+=local_gram_matrix[j*m_k+k];
270 int32_t* neighborhood_matrix =
SG_MALLOC(int32_t, N*k);
274 K_NEIGHBORHOOD_THREAD_PARAM* parameters =
SG_MALLOC(K_NEIGHBORHOOD_THREAD_PARAM, num_threads);
275 pthread_t* threads =
SG_MALLOC(pthread_t, num_threads);
277 pthread_attr_init(&attr);
278 pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
280 int32_t num_threads = 1;
282 CFibonacciHeap** heaps =
SG_MALLOC(CFibonacciHeap*, num_threads);
283 for (t=0; t<num_threads; t++)
284 heaps[t] =
new CFibonacciHeap(N);
287 for (t=0; t<num_threads; t++)
289 parameters[t].idx_start = t;
290 parameters[t].idx_step = num_threads;
291 parameters[t].idx_stop = N;
292 parameters[t].m_k = k;
294 parameters[t].heap = heaps[t];
295 parameters[t].neighborhood_matrix = neighborhood_matrix;
296 parameters[t].kernel_matrix = kernel_matrix.
matrix;
299 for (t=0; t<num_threads; t++)
300 pthread_join(threads[t], NULL);
301 pthread_attr_destroy(&attr);
305 K_NEIGHBORHOOD_THREAD_PARAM single_thread_param;
306 single_thread_param.idx_start = 0;
307 single_thread_param.idx_step = 1;
308 single_thread_param.idx_stop = N;
309 single_thread_param.m_k = k;
310 single_thread_param.N = N;
311 single_thread_param.heap = heaps[0]
312 single_thread_param.neighborhood_matrix = neighborhood_matrix;
313 single_thread_param.kernel_matrix = kernel_matrix.
matrix;
317 for (t=0; t<num_threads; t++)
326 K_NEIGHBORHOOD_THREAD_PARAM* parameters = (K_NEIGHBORHOOD_THREAD_PARAM*)p;
327 int32_t idx_start = parameters->idx_start;
328 int32_t idx_step = parameters->idx_step;
329 int32_t idx_stop = parameters->idx_stop;
330 int32_t N = parameters->N;
331 int32_t
m_k = parameters->m_k;
332 CFibonacciHeap* heap = parameters->heap;
333 const float64_t* kernel_matrix = parameters->kernel_matrix;
334 int32_t* neighborhood_matrix = parameters->neighborhood_matrix;
338 for (i=idx_start; i<idx_stop; i+=idx_step)
342 heap->insert(j,kernel_matrix[i*N+i]-2*kernel_matrix[i*N+j]+kernel_matrix[j*N+j]);
345 heap->extract_min(tmp);
347 for (j=0; j<
m_k; j++)
348 neighborhood_matrix[j*N+i] = heap->extract_min(tmp);