24 using namespace shogun;
26 #ifndef DOXYGEN_SHOULD_SKIP_THIS
27 struct HESSIANESTIMATION_THREAD_PARAM
40 const int32_t* neighborhood_matrix;
69 PTHREAD_LOCK_T* W_matrix_lock;
85 return "HessianLocallyLinearEmbedding";
95 SG_ERROR(
"K parameter should have value greater than 1+target dimensionality+dp.\n");
101 pthread_t* threads =
SG_MALLOC(pthread_t, num_threads);
102 HESSIANESTIMATION_THREAD_PARAM* parameters =
SG_MALLOC(HESSIANESTIMATION_THREAD_PARAM, num_threads);
105 int32_t num_threads = 1;
121 PTHREAD_LOCK_T W_matrix_lock;
123 PTHREAD_LOCK_INIT(&W_matrix_lock);
124 pthread_attr_init(&attr);
125 pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
127 for (t=0; t<num_threads; t++)
129 parameters[t].idx_start = t;
130 parameters[t].idx_step = num_threads;
131 parameters[t].idx_stop = N;
132 parameters[t].m_k =
m_k;
134 parameters[t].neighborhood_matrix = neighborhood_matrix.
matrix;
135 parameters[t].feature_matrix = feature_matrix.
matrix;
136 parameters[t].local_feature_matrix = local_feature_matrix + (
m_k*dim)*t;
138 parameters[t].mean_vector = mean_vector + dim*t;
139 parameters[t].s_values_vector = s_values_vector + dim*t;
140 parameters[t].tau = tau+tau_len*t;
141 parameters[t].tau_len = tau_len;
142 parameters[t].w_sum_vector = w_sum_vector + dp*t;
143 parameters[t].q_matrix = q_matrix + (
m_k*
m_k)*t;
144 parameters[t].W_matrix = W_matrix;
145 parameters[t].W_matrix_lock = &W_matrix_lock;
147 parameters[t].dim = dim;
148 parameters[t].actual_k = neighborhood_matrix.
num_rows;
151 for (t=0; t<num_threads; t++)
152 pthread_join(threads[t], NULL);
153 PTHREAD_LOCK_DESTROY(&W_matrix_lock);
157 HESSIANESTIMATION_THREAD_PARAM single_thread_param;
158 single_thread_param.idx_start = 0;
159 single_thread_param.idx_step = num_threads;
160 single_thread_param.idx_stop = N;
161 single_thread_param.m_k =
m_k;
162 single_thread_param.neighborhood_matrix = neighborhood_matrix.
matrix;
163 single_thread_param.feature_matrix = feature_matrix.
matrix;
164 single_thread_param.local_feature_matrix = local_feature_matrix;
165 single_thread_param.Yi_matrix = Yi_matrix;
166 single_thread_param.mean_vector = mean_vector;
167 single_thread_param.s_values_vector = s_values_vector;
168 single_thread_param.tau = tau;
169 single_thread_param.tau_len = tau_len;
170 single_thread_param.w_sum_vector = w_sum_vector;
171 single_thread_param.q_matrix = q_matrix;
172 single_thread_param.W_matrix = W_matrix;
173 single_thread_param.N = N;
174 single_thread_param.dim = dim;
175 single_thread_param.actual_k = neighborhood_matrix.
num_rows;
193 HESSIANESTIMATION_THREAD_PARAM* parameters = (HESSIANESTIMATION_THREAD_PARAM*)p;
194 int32_t idx_start = parameters->idx_start;
195 int32_t idx_step = parameters->idx_step;
196 int32_t idx_stop = parameters->idx_stop;
197 int32_t
m_k = parameters->m_k;
198 int32_t target_dim = parameters->target_dim;
199 const int32_t* neighborhood_matrix = parameters->neighborhood_matrix;
200 const float64_t* feature_matrix = parameters->feature_matrix;
201 float64_t* local_feature_matrix = parameters->local_feature_matrix;
202 float64_t* Yi_matrix = parameters->Yi_matrix;
203 float64_t* mean_vector = parameters->mean_vector;
204 float64_t* s_values_vector = parameters->s_values_vector;
206 int32_t tau_len = parameters->tau_len;
207 float64_t* w_sum_vector = parameters->w_sum_vector;
208 float64_t* q_matrix = parameters->q_matrix;
209 float64_t* W_matrix = parameters->W_matrix;
211 PTHREAD_LOCK_T* W_matrix_lock = parameters->W_matrix_lock;
215 int32_t dp = target_dim*(target_dim+1)/2;
216 int32_t actual_k = parameters->actual_k;
217 int32_t N = parameters->N;
218 int32_t dim = parameters->dim;
220 for (i=idx_start; i<idx_stop; i+=idx_step)
223 for (j=0; j<
m_k; j++)
227 memset(mean_vector,0,
sizeof(
float64_t)*dim);
230 for (j=0; j<
m_k; j++)
232 for (k=0; k<dim; k++)
233 local_feature_matrix[j*dim+k] = feature_matrix[neighborhood_matrix[i*actual_k+j]*dim+k];
235 cblas_daxpy(dim,1.0,local_feature_matrix+j*dim,1,mean_vector,1);
239 cblas_dscal(dim,1.0/m_k,mean_vector,1);
242 for (j=0; j<
m_k; j++)
243 cblas_daxpy(dim,-1.0,mean_vector,1,local_feature_matrix+j*dim,1);
247 wrap_dgesvd(
'N',
'O', dim,m_k,local_feature_matrix,dim,
249 NULL,1, NULL,1, &info);
253 for (j=0; j<target_dim; j++)
255 for (k=0; k<
m_k; k++)
256 Yi_matrix[(j+1)*m_k+k] = local_feature_matrix[k*dim+j];
262 for (j=0; j<target_dim; j++)
264 for (k=0; k<target_dim-j; k++)
266 for (l=0; l<
m_k; l++)
268 Yi_matrix[(ct+k+1+target_dim)*m_k+l] = Yi_matrix[(j+1)*m_k+l]*Yi_matrix[(j+k+1)*m_k+l];
271 ct += ct + target_dim - j;
275 wrap_dgeqrf(m_k,(1+target_dim+dp),Yi_matrix,m_k,tau,&info);
277 wrap_dorgqr(m_k,(1+target_dim+dp),tau_len,Yi_matrix,m_k,tau,&info);
280 float64_t* Pii = (Yi_matrix+m_k*(1+target_dim));
284 w_sum_vector[j] = 0.0;
285 for (k=0; k<
m_k; k++)
287 w_sum_vector[j] += Pii[j*m_k+k];
289 if (w_sum_vector[j]<0.001)
290 w_sum_vector[j] = 1.0;
291 for (k=0; k<
m_k; k++)
292 Pii[j*m_k+k] /= w_sum_vector[j];
295 cblas_dgemm(CblasColMajor,CblasNoTrans,CblasTrans,
301 PTHREAD_LOCK(W_matrix_lock);
303 for (j=0; j<
m_k; j++)
305 for (k=0; k<
m_k; k++)
306 W_matrix[N*neighborhood_matrix[i*actual_k+k]+neighborhood_matrix[i*actual_k+j]] += q_matrix[j*m_k+k];
309 PTHREAD_UNLOCK(W_matrix_lock);