Mercurial > hg > smallbox
comparison util/ksvd utils/ompbox utils/ompcoreGabor.c @ 137:9207d56c5547 ivand_dev
New ompbox in utils for testing purposes
author | Ivan Damnjanovic lnx <ivan.damnjanovic@eecs.qmul.ac.uk> |
---|---|
date | Thu, 21 Jul 2011 14:07:41 +0100 |
parents | |
children | 4bd6856a7128 |
comparison
equal
deleted
inserted
replaced
136:1334d2302dd9 | 137:9207d56c5547 |
---|---|
1 /************************************************************************** | |
2 * | |
3 * File name: ompcoreGabor.c | |
4 * | |
5 * Ron Rubinstein | |
6 * Computer Science Department | |
7 * Technion, Haifa 32000 Israel | |
8 * ronrubin@cs | |
9 * | |
10 * Last Updated: 25.8.2009 | |
11 * | |
12 * Modified by Ivan damnjanovic July 2011 | |
13 * Takes to atoms per iteration. It should be used for Gabor dictionaries | |
14 * as specified in | |
15 * "Audio Inpainting" Amir Adler, Valentin Emiya, Maria G. Jafari, | |
16 * Michael Elad, Remi Gribonval and Mark D. Plumbley | |
17 * Draft version: March 6, 2011 | |
18 * | |
19 *************************************************************************/ | |
20 | |
21 | |
22 #include "ompcoreGabor.h" | |
23 #include "omputils.h" | |
24 #include "ompprof.h" | |
25 #include "myblas.h" | |
26 #include <math.h> | |
27 #include <string.h> | |
28 | |
29 | |
30 | |
31 /****************************************************************************** | |
32 * * | |
33 * Batch-OMP Implementation * | |
34 * * | |
35 ******************************************************************************/ | |
36 | |
37 mxArray* ompcoreGabor(double D[], double x[], double DtX[], double XtX[], double G[], mwSize n, mwSize m, mwSize L, | |
38 int T, double eps, int gamma_mode, int profile, double msg_delta, int erroromp) | |
39 { | |
40 | |
41 profdata pd; | |
42 mxArray *Gamma; | |
43 mwIndex i, j, k, signum, pos, *ind, *gammaIr, *gammaJc, gamma_count; | |
44 mwSize allocated_coefs, allocated_cols; | |
45 int DtX_specified, XtX_specified, batchomp, standardomp, *selected_atoms; | |
46 double *proj, *proj1, *proj2, *D1, *D2, *D1D2, *n12, *alpha, *beta, *error; | |
47 double *r, *Lchol, *c, *Gsub, *Dsub, sum, *gammaPr, *tempvec1, *tempvec2; | |
48 double eps2, resnorm, delta, deltaprev, secs_remain; | |
49 int mins_remain, hrs_remain; | |
50 clock_t lastprint_time, starttime; | |
51 | |
52 | |
53 | |
54 /*** status flags ***/ | |
55 | |
56 DtX_specified = (DtX!=0); /* indicates whether D'*x was provided */ | |
57 XtX_specified = (XtX!=0); /* indicates whether sum(x.*x) was provided */ | |
58 | |
59 standardomp = (G==0); /* batch-omp or standard omp are selected depending on availability of G */ | |
60 batchomp = !standardomp; | |
61 | |
62 | |
63 | |
64 /*** allocate output matrix ***/ | |
65 | |
66 | |
67 if (gamma_mode == FULL_GAMMA) { | |
68 | |
69 /* allocate full matrix of size m X L */ | |
70 | |
71 Gamma = mxCreateDoubleMatrix(m, L, mxREAL); | |
72 gammaPr = mxGetPr(Gamma); | |
73 gammaIr = 0; | |
74 gammaJc = 0; | |
75 } | |
76 else { | |
77 | |
78 /* allocate sparse matrix with room for allocated_coefs nonzeros */ | |
79 | |
80 /* for error-omp, begin with L*sqrt(n)/2 allocated nonzeros, otherwise allocate L*T nonzeros */ | |
81 allocated_coefs = erroromp ? (mwSize)(ceil(L*sqrt((double)n)/2.0) + 1.01) : L*T; | |
82 Gamma = mxCreateSparse(m, L, allocated_coefs, mxREAL); | |
83 gammaPr = mxGetPr(Gamma); | |
84 gammaIr = mxGetIr(Gamma); | |
85 gammaJc = mxGetJc(Gamma); | |
86 gamma_count = 0; | |
87 gammaJc[0] = 0; | |
88 } | |
89 | |
90 | |
91 /*** helper arrays ***/ | |
92 /* Ivan Damnjanovic July 2011*/ | |
93 proj = (double*)mxMalloc(m*sizeof(double)); | |
94 proj1 = (double*)mxMalloc(m/2*sizeof(double)); | |
95 proj2 = (double*)mxMalloc(m/2*sizeof(double)); | |
96 D1 = (double*)mxMalloc(n*m/2*sizeof(double)); | |
97 D2 = (double*)mxMalloc(n*m/2*sizeof(double)); | |
98 memcpy(D1, D , n*m/2*sizeof(double)); | |
99 memcpy(D2, D+n*m/2, n*m/2*sizeof(double)); | |
100 | |
101 D1D2 = (double*)mxMalloc(m/2*sizeof(double)); | |
102 n12 = (double*)mxMalloc(m/2*sizeof(double)); | |
103 | |
104 vec_smult(1,D2, D1, n*m/2); | |
105 | |
106 for (i=0; i<m/2; i++) { | |
107 for (j=0; j<n; j++) { | |
108 D1D2[i] += D1[i*n+j]; | |
109 } | |
110 n12[i]=1/(1-D1D2[i]); | |
111 } | |
112 | |
113 memcpy(D1, D , n*m/2*sizeof(double)); | |
114 | |
115 alpha = (double*)mxMalloc(m/2*sizeof(double)); /* contains D'*residual */ | |
116 beta = (double*)mxMalloc(m/2*sizeof(double)); | |
117 error = (double*)mxMalloc(m/2*sizeof(double)); | |
118 | |
119 ind = (mwIndex*)mxMalloc(m*sizeof(mwIndex)); /* indices of selected atoms */ | |
120 selected_atoms = (int*)mxMalloc(m*sizeof(int)); /* binary array with 1's for selected atoms */ | |
121 c = (double*)mxMalloc(n*sizeof(double)); /* orthogonal projection result */ | |
122 | |
123 /* current number of columns in Dsub / Gsub / Lchol */ | |
124 allocated_cols = erroromp ? (mwSize)(ceil(sqrt((double)n)/2.0) + 1.01) : T; | |
125 | |
126 /* Cholesky decomposition of D_I'*D_I */ | |
127 Lchol = (double*)mxMalloc(n*allocated_cols*sizeof(double)); | |
128 | |
129 /* temporary vectors for various computations */ | |
130 tempvec1 = (double*)mxMalloc(m*sizeof(double)); | |
131 tempvec2 = (double*)mxMalloc(m*sizeof(double)); | |
132 | |
133 if (batchomp) { | |
134 /* matrix containing G(:,ind) - the columns of G corresponding to the selected atoms, in order of selection */ | |
135 Gsub = (double*)mxMalloc(m*allocated_cols*sizeof(double)); | |
136 } | |
137 else { | |
138 /* matrix containing D(:,ind) - the selected atoms from D, in order of selection */ | |
139 Dsub = (double*)mxMalloc(n*allocated_cols*sizeof(double)); | |
140 | |
141 /* stores the residual */ | |
142 r = (double*)mxMalloc(n*sizeof(double)); | |
143 } | |
144 | |
145 if (!DtX_specified) { | |
146 /* contains D'*x for the current signal */ | |
147 DtX = (double*)mxMalloc(m*sizeof(double)); | |
148 } | |
149 | |
150 | |
151 | |
152 /*** initializations for error omp ***/ | |
153 | |
154 if (erroromp) { | |
155 eps2 = eps*eps; /* compute eps^2 */ | |
156 if (T<0 || T>n) { /* unspecified max atom num - set max atoms to n */ | |
157 T = n; | |
158 } | |
159 } | |
160 | |
161 | |
162 | |
163 /*** initialize timers ***/ | |
164 | |
165 initprofdata(&pd); /* initialize profiling counters */ | |
166 starttime = clock(); /* record starting time for eta computations */ | |
167 lastprint_time = starttime; /* time of last status display */ | |
168 | |
169 | |
170 | |
171 /********************** perform omp for each signal **********************/ | |
172 | |
173 | |
174 | |
175 for (signum=0; signum<L; ++signum) { | |
176 | |
177 | |
178 /* initialize residual norm and deltaprev for error-omp */ | |
179 | |
180 if (erroromp) { | |
181 if (XtX_specified) { | |
182 resnorm = XtX[signum]; | |
183 } | |
184 else { | |
185 resnorm = dotprod(x+n*signum, x+n*signum, n); | |
186 addproftime(&pd, XtX_TIME); | |
187 } | |
188 deltaprev = 0; /* delta tracks the value of gamma'*G*gamma */ | |
189 } | |
190 else { | |
191 /* ignore residual norm stopping criterion */ | |
192 eps2 = 0; | |
193 resnorm = 1; | |
194 } | |
195 | |
196 | |
197 if (resnorm>eps2 && T>0) { | |
198 | |
199 /* compute DtX */ | |
200 | |
201 if (!DtX_specified) { | |
202 matT_vec(1, D, x+n*signum, DtX, n, m); | |
203 addproftime(&pd, DtX_TIME); | |
204 memcpy(r , x+n*signum, n*sizeof(double)); | |
205 } | |
206 | |
207 | |
208 /* initialize projections to D1 and D2 := DtX */ | |
209 | |
210 memcpy(proj, DtX + m*signum*DtX_specified, m*sizeof(double)); | |
211 | |
212 | |
213 /* mark all atoms as unselected */ | |
214 | |
215 for (i=0; i<m; ++i) { | |
216 selected_atoms[i] = 0; | |
217 } | |
218 | |
219 } | |
220 | |
221 | |
222 /* main loop */ | |
223 | |
224 i=0; | |
225 while (resnorm>eps2 && i<T) { | |
226 | |
227 /* index of next atom */ | |
228 memcpy(proj1, proj, m/2*sizeof(double)); | |
229 memcpy(proj2, proj + m/2, m/2*sizeof(double)); | |
230 for (k=0; k<m/2; k++){ | |
231 alpha[k] = (proj1[k] - D1D2[k]*proj2[k])*n12[k]; | |
232 beta[k] = (proj2[k] - D1D2[k]*proj1[k])*n12[k]; | |
233 } | |
234 for (k=0; k<m/2; k++){ | |
235 for (j=0; j<n; j++){ | |
236 error[k]+= (abs(r[j])-D1[k*n+j]*alpha[k]-D2[(k+m/2)*n+j]*beta[k])*(abs(r[j])-D1[k*n+j]*alpha[k]-D2[(k+m/2)*n+j]*beta[k]); | |
237 } | |
238 } | |
239 pos = maxabs(error, m/2); | |
240 addproftime(&pd, MAXABS_TIME); | |
241 | |
242 | |
243 /* stop criterion: selected same atom twice, or inner product too small */ | |
244 | |
245 if (selected_atoms[pos] || alpha[pos]*alpha[pos]<1e-14) { | |
246 break; | |
247 } | |
248 | |
249 for (k=0;k<2;k++){ | |
250 /* mark selected atom */ | |
251 | |
252 ind[i] = pos+k*m/2; | |
253 selected_atoms[pos+k*m/2] = 1; | |
254 | |
255 | |
256 /* matrix reallocation */ | |
257 | |
258 if (erroromp && i>=allocated_cols) { | |
259 | |
260 allocated_cols = (mwSize)(ceil(allocated_cols*MAT_INC_FACTOR) + 1.01); | |
261 | |
262 Lchol = (double*)mxRealloc(Lchol,n*allocated_cols*sizeof(double)); | |
263 | |
264 batchomp ? (Gsub = (double*)mxRealloc(Gsub,m*allocated_cols*sizeof(double))) : | |
265 (Dsub = (double*)mxRealloc(Dsub,n*allocated_cols*sizeof(double))) ; | |
266 } | |
267 | |
268 | |
269 /* append column to Gsub or Dsub */ | |
270 | |
271 if (batchomp) { | |
272 memcpy(Gsub+i*m, G+(pos+k*m/2)*m, m*sizeof(double)); | |
273 } | |
274 else { | |
275 memcpy(Dsub+(i)*n, D+(pos+k*m/2)*n, n*sizeof(double)); | |
276 } | |
277 | |
278 | |
279 /*** Cholesky update ***/ | |
280 | |
281 if (i==0) { | |
282 *Lchol = 1; | |
283 } | |
284 else { | |
285 | |
286 /* incremental Cholesky decomposition: compute next row of Lchol */ | |
287 | |
288 if (standardomp) { | |
289 matT_vec(1, Dsub, D+n*(pos+k*m/2), tempvec1, n, i); /* compute tempvec1 := Dsub'*d where d is new atom */ | |
290 addproftime(&pd, DtD_TIME); | |
291 } | |
292 else { | |
293 vec_assign(tempvec1, Gsub+i*m, ind, i); /* extract tempvec1 := Gsub(ind,i) */ | |
294 } | |
295 backsubst('L', Lchol, tempvec1, tempvec2, n, i); /* compute tempvec2 = Lchol \ tempvec1 */ | |
296 for (j=0; j<i; ++j) { /* write tempvec2 to end of Lchol */ | |
297 Lchol[j*n+i] = tempvec2[j]; | |
298 } | |
299 | |
300 /* compute Lchol(i,i) */ | |
301 sum = 0; | |
302 for (j=0; j<i; ++j) { /* compute sum of squares of last row without Lchol(i,i) */ | |
303 sum += SQR(Lchol[j*n+i]); | |
304 } | |
305 if ( (1-sum) <= 1e-14 ) { /* Lchol(i,i) is zero => selected atoms are dependent */ | |
306 break; | |
307 } | |
308 Lchol[i*n+i] = sqrt(1-sum); | |
309 } | |
310 | |
311 addproftime(&pd, LCHOL_TIME); | |
312 | |
313 i++; | |
314 | |
315 } | |
316 /* perform orthogonal projection and compute sparse coefficients */ | |
317 | |
318 vec_assign(tempvec1, DtX + m*signum*DtX_specified, ind, i); /* extract tempvec1 = DtX(ind) */ | |
319 cholsolve('L', Lchol, tempvec1, c, n, i); /* solve LL'c = tempvec1 for c */ | |
320 addproftime(&pd, COMPCOEF_TIME); | |
321 | |
322 | |
323 /* update alpha = D'*residual */ | |
324 | |
325 if (standardomp) { | |
326 mat_vec(-1, Dsub, c, r, n, i); /* compute r := -Dsub*c */ | |
327 vec_sum(1, x+n*signum, r, n); /* compute r := x+r */ | |
328 | |
329 | |
330 /*memcpy(r, x+n*signum, n*sizeof(double)); /* assign r := x */ | |
331 /*mat_vec1(-1, Dsub, c, 1, r, n, i); /* compute r := r-Dsub*c */ | |
332 | |
333 addproftime(&pd, COMPRES_TIME); | |
334 matT_vec(1, D, r, proj, n, m); /* compute proj := D'*r */ | |
335 addproftime(&pd, DtR_TIME); | |
336 | |
337 /* update residual norm */ | |
338 if (erroromp) { | |
339 resnorm = dotprod(r, r, n); | |
340 addproftime(&pd, UPDATE_RESNORM_TIME); | |
341 } | |
342 } | |
343 else { | |
344 mat_vec(1, Gsub, c, tempvec1, m, i); /* compute tempvec1 := Gsub*c */ | |
345 memcpy(proj, DtX + m*signum*DtX_specified, m*sizeof(double)); /* set proj = D'*x */ | |
346 vec_sum(-1, tempvec1, proj, m); /* compute proj := proj - tempvec1 */ | |
347 addproftime(&pd, UPDATE_DtR_TIME); | |
348 | |
349 /* update residual norm */ | |
350 if (erroromp) { | |
351 vec_assign(tempvec2, tempvec1, ind, i); /* assign tempvec2 := tempvec1(ind) */ | |
352 delta = dotprod(c,tempvec2,i); /* compute c'*tempvec2 */ | |
353 resnorm = resnorm - delta + deltaprev; /* residual norm update */ | |
354 deltaprev = delta; | |
355 addproftime(&pd, UPDATE_RESNORM_TIME); | |
356 } | |
357 } | |
358 } | |
359 | |
360 | |
361 /*** generate output vector gamma ***/ | |
362 | |
363 if (gamma_mode == FULL_GAMMA) { /* write the coefs in c to their correct positions in gamma */ | |
364 for (j=0; j<i; ++j) { | |
365 gammaPr[m*signum + ind[j]] = c[j]; | |
366 } | |
367 } | |
368 else { | |
369 /* sort the coefs by index before writing them to gamma */ | |
370 quicksort(ind,c,i); | |
371 addproftime(&pd, INDEXSORT_TIME); | |
372 | |
373 /* gamma is full - reallocate */ | |
374 if (gamma_count+i >= allocated_coefs) { | |
375 | |
376 while(gamma_count+i >= allocated_coefs) { | |
377 allocated_coefs = (mwSize)(ceil(GAMMA_INC_FACTOR*allocated_coefs) + 1.01); | |
378 } | |
379 | |
380 mxSetNzmax(Gamma, allocated_coefs); | |
381 mxSetPr(Gamma, mxRealloc(gammaPr, allocated_coefs*sizeof(double))); | |
382 mxSetIr(Gamma, mxRealloc(gammaIr, allocated_coefs*sizeof(mwIndex))); | |
383 | |
384 gammaPr = mxGetPr(Gamma); | |
385 gammaIr = mxGetIr(Gamma); | |
386 } | |
387 | |
388 /* append coefs to gamma and update the indices */ | |
389 for (j=0; j<i; ++j) { | |
390 gammaPr[gamma_count] = c[j]; | |
391 gammaIr[gamma_count] = ind[j]; | |
392 gamma_count++; | |
393 } | |
394 gammaJc[signum+1] = gammaJc[signum] + i; | |
395 } | |
396 | |
397 | |
398 | |
399 /*** display status messages ***/ | |
400 | |
401 if (msg_delta>0 && (clock()-lastprint_time)/(double)CLOCKS_PER_SEC >= msg_delta) | |
402 { | |
403 lastprint_time = clock(); | |
404 | |
405 /* estimated remainig time */ | |
406 secs2hms( ((L-signum-1)/(double)(signum+1)) * ((lastprint_time-starttime)/(double)CLOCKS_PER_SEC) , | |
407 &hrs_remain, &mins_remain, &secs_remain); | |
408 | |
409 mexPrintf("omp: signal %d / %d, estimated remaining time: %02d:%02d:%05.2f\n", | |
410 signum+1, L, hrs_remain, mins_remain, secs_remain); | |
411 mexEvalString("drawnow;"); | |
412 } | |
413 | |
414 } | |
415 | |
416 /* end omp */ | |
417 | |
418 | |
419 | |
420 /*** print final messages ***/ | |
421 | |
422 if (msg_delta>0) { | |
423 mexPrintf("omp: signal %d / %d\n", signum, L); | |
424 } | |
425 | |
426 if (profile) { | |
427 printprofinfo(&pd, erroromp, batchomp, L); | |
428 } | |
429 | |
430 | |
431 | |
432 /* free memory */ | |
433 | |
434 if (!DtX_specified) { | |
435 mxFree(DtX); | |
436 } | |
437 if (standardomp) { | |
438 mxFree(r); | |
439 mxFree(Dsub); | |
440 } | |
441 else { | |
442 mxFree(Gsub); | |
443 } | |
444 mxFree(tempvec2); | |
445 mxFree(tempvec1); | |
446 mxFree(Lchol); | |
447 mxFree(c); | |
448 mxFree(selected_atoms); | |
449 mxFree(ind); | |
450 mxFree(proj); | |
451 mxFree(proj1); | |
452 mxFree(proj2); | |
453 mxFree(D1); | |
454 mxFree(D2); | |
455 mxFree(D1D2); | |
456 mxFree(n12); | |
457 mxFree(alpha); | |
458 mxFree(beta); | |
459 mxFree(error); | |
460 | |
461 return Gamma; | |
462 } |