Mercurial > hg > smallbox
view solvers/SMALL_ompGabor/ompcoreGabor.c @ 236:5f4e47b78f2b ver_2.0_alpha1
Added example.
author | luisf <luis.figueira@eecs.qmul.ac.uk> |
---|---|
date | Thu, 19 Apr 2012 17:59:08 +0100 |
parents | 31d2864dfdd4 |
children |
line wrap: on
line source
/************************************************************************** * * File name: ompcoreGabor.c * * Ron Rubinstein * Computer Science Department * Technion, Haifa 32000 Israel * ronrubin@cs * * Last Updated: 25.8.2009 * * Modified by Ivan damnjanovic July 2011 * Takes to atoms per iteration. It should be used for Gabor dictionaries * as specified in * "Audio Inpainting" Amir Adler, Valentin Emiya, Maria G. Jafari, * Michael Elad, Remi Gribonval and Mark D. Plumbley * Draft version: March 6, 2011 * *************************************************************************/ #include "ompcoreGabor.h" #include "omputils.h" #include "ompprof.h" #include "myblas.h" #include <math.h> #include <string.h> /****************************************************************************** * * * Batch-OMP Implementation * * * ******************************************************************************/ mxArray* ompcoreGabor(double D[], double x[], double DtX[], double XtX[], double G[], mwSize n, mwSize m, mwSize L, int T, double eps, int gamma_mode, int profile, double msg_delta, int erroromp) { profdata pd; mxArray *Gamma; mwIndex i, j, k, signum, pos, *ind, *gammaIr, *gammaJc, gamma_count; mwSize allocated_coefs, allocated_cols; int DtX_specified, XtX_specified, batchomp, standardomp, *selected_atoms; double *proj, *proj1, *proj2, *D1, *D2, *D1D2, *n12, *alpha, *beta, *error; double *r, *Lchol, *c, *Gsub, *Dsub, sum, *gammaPr, *tempvec1, *tempvec2; double eps2, resnorm, delta, deltaprev, secs_remain; int mins_remain, hrs_remain; clock_t lastprint_time, starttime; /*** status flags ***/ DtX_specified = (DtX!=0); /* indicates whether D'*x was provided */ XtX_specified = (XtX!=0); /* indicates whether sum(x.*x) was provided */ standardomp = (G==0); /* batch-omp or standard omp are selected depending on availability of G */ batchomp = !standardomp; /*** allocate output matrix ***/ if (gamma_mode == FULL_GAMMA) { /* allocate full matrix of size m X L */ Gamma = mxCreateDoubleMatrix(m, L, mxREAL); gammaPr = mxGetPr(Gamma); gammaIr = 0; gammaJc = 0; } else { /* allocate sparse matrix with room for allocated_coefs nonzeros */ /* for error-omp, begin with L*sqrt(n)/2 allocated nonzeros, otherwise allocate L*T nonzeros */ allocated_coefs = erroromp ? (mwSize)(ceil(L*sqrt((double)n)/2.0) + 1.01) : L*T; Gamma = mxCreateSparse(m, L, allocated_coefs, mxREAL); gammaPr = mxGetPr(Gamma); gammaIr = mxGetIr(Gamma); gammaJc = mxGetJc(Gamma); gamma_count = 0; gammaJc[0] = 0; } /*** helper arrays ***/ /* Ivan Damnjanovic July 2011*/ proj = (double*)mxMalloc(m*sizeof(double)); proj1 = (double*)mxMalloc(m/2*sizeof(double)); proj2 = (double*)mxMalloc(m/2*sizeof(double)); D1 = (double*)mxMalloc(n*m/2*sizeof(double)); D2 = (double*)mxMalloc(n*m/2*sizeof(double)); memcpy(D1, D , n*m/2*sizeof(double)); memcpy(D2, D+n*m/2, n*m/2*sizeof(double)); D1D2 = (double*)mxMalloc(m/2*sizeof(double)); n12 = (double*)mxMalloc(m/2*sizeof(double)); vec_smult(1,D2, D1, n*m/2); for (i=0; i<m/2; i++) { D1D2[i]=0; n12[i]=0; for (j=0; j<n; j++) { D1D2[i] += D1[i*n+j]; } n12[i]=1/(1-D1D2[i]*D1D2[i]); } memcpy(D1, D , n*m/2*sizeof(double)); alpha = (double*)mxMalloc(m/2*sizeof(double)); /* contains D'*residual */ beta = (double*)mxMalloc(m/2*sizeof(double)); error = (double*)mxMalloc(m/2*sizeof(double)); ind = (mwIndex*)mxMalloc(m*sizeof(mwIndex)); /* indices of selected atoms */ selected_atoms = (int*)mxMalloc(m*sizeof(int)); /* binary array with 1's for selected atoms */ c = (double*)mxMalloc(n*sizeof(double)); /* orthogonal projection result */ /* current number of columns in Dsub / Gsub / Lchol */ allocated_cols = erroromp ? (mwSize)(ceil(sqrt((double)n)/2.0) + 1.01) : T; /* Cholesky decomposition of D_I'*D_I */ Lchol = (double*)mxMalloc(n*allocated_cols*sizeof(double)); /* temporary vectors for various computations */ tempvec1 = (double*)mxMalloc(m*sizeof(double)); tempvec2 = (double*)mxMalloc(m*sizeof(double)); if (batchomp) { /* matrix containing G(:,ind) - the columns of G corresponding to the selected atoms, in order of selection */ Gsub = (double*)mxMalloc(m*allocated_cols*sizeof(double)); } else { /* matrix containing D(:,ind) - the selected atoms from D, in order of selection */ Dsub = (double*)mxMalloc(n*allocated_cols*sizeof(double)); /* stores the residual */ r = (double*)mxMalloc(n*sizeof(double)); } if (!DtX_specified) { /* contains D'*x for the current signal */ DtX = (double*)mxMalloc(m*sizeof(double)); } /*** initializations for error omp ***/ if (erroromp) { eps2 = eps*eps; /* compute eps^2 */ if (T<0 || T>n) { /* unspecified max atom num - set max atoms to n */ T = n; } } /*** initialize timers ***/ initprofdata(&pd); /* initialize profiling counters */ starttime = clock(); /* record starting time for eta computations */ lastprint_time = starttime; /* time of last status display */ /********************** perform omp for each signal **********************/ for (signum=0; signum<L; ++signum) { /* initialize residual norm and deltaprev for error-omp */ if (erroromp) { if (XtX_specified) { resnorm = XtX[signum]; } else { resnorm = dotprod(x+n*signum, x+n*signum, n); addproftime(&pd, XtX_TIME); } deltaprev = 0; /* delta tracks the value of gamma'*G*gamma */ } else { /* ignore residual norm stopping criterion */ eps2 = 0; resnorm = 1; } if (resnorm>eps2 && T>0) { /* compute DtX */ if (!DtX_specified) { matT_vec(1, D, x+n*signum, DtX, n, m); addproftime(&pd, DtX_TIME); memcpy(r , x+n*signum, n*sizeof(double)); } /* initialize projections to D1 and D2 := DtX */ memcpy(proj, DtX + m*signum*DtX_specified, m*sizeof(double)); /* mark all atoms as unselected */ for (i=0; i<m; ++i) { selected_atoms[i] = 0; } } /* main loop */ i=0; while (resnorm>eps2 && i<T) { /* index of next atom */ memcpy(proj1, proj, m/2*sizeof(double)); memcpy(proj2, proj + m/2, m/2*sizeof(double)); for (k=0; k<m/2; k++){ alpha[k] = (proj1[k] - D1D2[k]*proj2[k])*n12[k]; beta[k] = (proj2[k] - D1D2[k]*proj1[k])*n12[k]; } for (k=0; k<m/2; k++){ error[k]=0; for (j=0; j<n; j++){ error[k]+= (abs(r[j])-D1[k*n+j]*alpha[k]-D2[k*n+j]*beta[k])*(abs(r[j])-D1[k*n+j]*alpha[k]-D2[k*n+j]*beta[k]); } } pos = maxabs(error, m/2); addproftime(&pd, MAXABS_TIME); /* stop criterion: selected same atom twice, or inner product too small */ if (selected_atoms[pos] || alpha[pos]*alpha[pos]<1e-14) { break; } for (k=0;k<2;k++){ /* mark selected atom */ ind[i] = pos+k*m/2; selected_atoms[pos+k*m/2] = 1; /* matrix reallocation */ if (erroromp && i>=allocated_cols) { allocated_cols = (mwSize)(ceil(allocated_cols*MAT_INC_FACTOR) + 1.01); Lchol = (double*)mxRealloc(Lchol,n*allocated_cols*sizeof(double)); batchomp ? (Gsub = (double*)mxRealloc(Gsub,m*allocated_cols*sizeof(double))) : (Dsub = (double*)mxRealloc(Dsub,n*allocated_cols*sizeof(double))) ; } /* append column to Gsub or Dsub */ if (batchomp) { memcpy(Gsub+i*m, G+(pos+k*m/2)*m, m*sizeof(double)); } else { memcpy(Dsub+(i)*n, D+(pos+k*m/2)*n, n*sizeof(double)); } /*** Cholesky update ***/ if (i==0) { *Lchol = 1; } else { /* incremental Cholesky decomposition: compute next row of Lchol */ if (standardomp) { matT_vec(1, Dsub, D+n*(pos+k*m/2), tempvec1, n, i); /* compute tempvec1 := Dsub'*d where d is new atom */ addproftime(&pd, DtD_TIME); } else { vec_assign(tempvec1, Gsub+i*m, ind, i); /* extract tempvec1 := Gsub(ind,i) */ } backsubst('L', Lchol, tempvec1, tempvec2, n, i); /* compute tempvec2 = Lchol \ tempvec1 */ for (j=0; j<i; ++j) { /* write tempvec2 to end of Lchol */ Lchol[j*n+i] = tempvec2[j]; } /* compute Lchol(i,i) */ sum = 0; for (j=0; j<i; ++j) { /* compute sum of squares of last row without Lchol(i,i) */ sum += SQR(Lchol[j*n+i]); } if ( (1-sum) <= 1e-14 ) { /* Lchol(i,i) is zero => selected atoms are dependent */ break; } Lchol[i*n+i] = sqrt(1-sum); } addproftime(&pd, LCHOL_TIME); i++; } /* perform orthogonal projection and compute sparse coefficients */ vec_assign(tempvec1, DtX + m*signum*DtX_specified, ind, i); /* extract tempvec1 = DtX(ind) */ cholsolve('L', Lchol, tempvec1, c, n, i); /* solve LL'c = tempvec1 for c */ addproftime(&pd, COMPCOEF_TIME); /* update alpha = D'*residual */ if (standardomp) { mat_vec(-1, Dsub, c, r, n, i); /* compute r := -Dsub*c */ vec_sum(1, x+n*signum, r, n); /* compute r := x+r */ /*memcpy(r, x+n*signum, n*sizeof(double)); /* assign r := x */ /*mat_vec1(-1, Dsub, c, 1, r, n, i); /* compute r := r-Dsub*c */ addproftime(&pd, COMPRES_TIME); matT_vec(1, D, r, proj, n, m); /* compute proj := D'*r */ addproftime(&pd, DtR_TIME); /* update residual norm */ if (erroromp) { resnorm = dotprod(r, r, n); addproftime(&pd, UPDATE_RESNORM_TIME); } } else { mat_vec(1, Gsub, c, tempvec1, m, i); /* compute tempvec1 := Gsub*c */ memcpy(proj, DtX + m*signum*DtX_specified, m*sizeof(double)); /* set proj = D'*x */ vec_sum(-1, tempvec1, proj, m); /* compute proj := proj - tempvec1 */ addproftime(&pd, UPDATE_DtR_TIME); /* update residual norm */ if (erroromp) { vec_assign(tempvec2, tempvec1, ind, i); /* assign tempvec2 := tempvec1(ind) */ delta = dotprod(c,tempvec2,i); /* compute c'*tempvec2 */ resnorm = resnorm - delta + deltaprev; /* residual norm update */ deltaprev = delta; addproftime(&pd, UPDATE_RESNORM_TIME); } } } /*** generate output vector gamma ***/ if (gamma_mode == FULL_GAMMA) { /* write the coefs in c to their correct positions in gamma */ for (j=0; j<i; ++j) { gammaPr[m*signum + ind[j]] = c[j]; } } else { /* sort the coefs by index before writing them to gamma */ quicksort(ind,c,i); addproftime(&pd, INDEXSORT_TIME); /* gamma is full - reallocate */ if (gamma_count+i >= allocated_coefs) { while(gamma_count+i >= allocated_coefs) { allocated_coefs = (mwSize)(ceil(GAMMA_INC_FACTOR*allocated_coefs) + 1.01); } mxSetNzmax(Gamma, allocated_coefs); mxSetPr(Gamma, mxRealloc(gammaPr, allocated_coefs*sizeof(double))); mxSetIr(Gamma, mxRealloc(gammaIr, allocated_coefs*sizeof(mwIndex))); gammaPr = mxGetPr(Gamma); gammaIr = mxGetIr(Gamma); } /* append coefs to gamma and update the indices */ for (j=0; j<i; ++j) { gammaPr[gamma_count] = c[j]; gammaIr[gamma_count] = ind[j]; gamma_count++; } gammaJc[signum+1] = gammaJc[signum] + i; } /*** display status messages ***/ if (msg_delta>0 && (clock()-lastprint_time)/(double)CLOCKS_PER_SEC >= msg_delta) { lastprint_time = clock(); /* estimated remainig time */ secs2hms( ((L-signum-1)/(double)(signum+1)) * ((lastprint_time-starttime)/(double)CLOCKS_PER_SEC) , &hrs_remain, &mins_remain, &secs_remain); mexPrintf("omp: signal %d / %d, estimated remaining time: %02d:%02d:%05.2f\n", signum+1, L, hrs_remain, mins_remain, secs_remain); mexEvalString("drawnow;"); } } /* end omp */ /*** print final messages ***/ if (msg_delta>0) { mexPrintf("omp: signal %d / %d\n", signum, L); } if (profile) { printprofinfo(&pd, erroromp, batchomp, L); } /* free memory */ if (!DtX_specified) { mxFree(DtX); } if (standardomp) { mxFree(r); mxFree(Dsub); } else { mxFree(Gsub); } mxFree(tempvec2); mxFree(tempvec1); mxFree(Lchol); mxFree(c); mxFree(selected_atoms); mxFree(ind); mxFree(proj); mxFree(proj1); mxFree(proj2); mxFree(D1); mxFree(D2); mxFree(D1D2); mxFree(n12); mxFree(alpha); mxFree(beta); mxFree(error); return Gamma; }