changeset 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 1334d2302dd9
children 56d719a5fd31
files Problems/generateAudioDeclippingProblem.m examples/AudioInpainting/Audio_Declipping_Example.m util/ksvd utils/ompbox utils/make.m util/ksvd utils/ompbox utils/mexutils.c util/ksvd utils/ompbox utils/mexutils.h util/ksvd utils/ompbox utils/myblas.c util/ksvd utils/ompbox utils/myblas.h util/ksvd utils/ompbox utils/omp2Gabor.m util/ksvd utils/ompbox utils/omp2mex.c util/ksvd utils/ompbox utils/omp2mex.m util/ksvd utils/ompbox utils/omp2mexGabor.c util/ksvd utils/ompbox utils/omp2mexGabor.m util/ksvd utils/ompbox utils/ompGabor.m util/ksvd utils/ompbox utils/ompcore.c util/ksvd utils/ompbox utils/ompcore.h util/ksvd utils/ompbox utils/ompcoreGabor.c util/ksvd utils/ompbox utils/ompcoreGabor.h util/ksvd utils/ompbox utils/ompmex.c util/ksvd utils/ompbox utils/ompmex.m util/ksvd utils/ompbox utils/ompmexGabor.c util/ksvd utils/ompbox utils/ompmexGabor.m util/ksvd utils/ompbox utils/ompprof.c util/ksvd utils/ompbox utils/ompprof.h util/ksvd utils/ompbox utils/omputils.c util/ksvd utils/ompbox utils/omputils.h util/ksvd utils/ompbox utils/printf.m
diffstat 26 files changed, 3944 insertions(+), 13 deletions(-) [+]
line wrap: on
line diff
--- a/Problems/generateAudioDeclippingProblem.m	Thu Jul 14 16:26:07 2011 +0100
+++ b/Problems/generateAudioDeclippingProblem.m	Thu Jul 21 14:07:41 2011 +0100
@@ -88,7 +88,7 @@
 
 x_clip = im2colstep(problemData.x,[windowSize 1],[overlap*windowSize 1]);
 x_clip= diag(wa(windowSize)) * x_clip;
-blkMask=im2colstep(double(~problemData.IMiss),[256 1],[128 1]);
+blkMask=im2colstep(double(~problemData.IMiss),[windowSize 1],[overlap*windowSize 1]);
 
 %% Building dictionary
 if ~exist( 'redundancyFactor', 'var' ) || isempty(redundancyFactor), redundancyFactor = 2; end % Weighting window for dictionary atoms
--- a/examples/AudioInpainting/Audio_Declipping_Example.m	Thu Jul 14 16:26:07 2011 +0100
+++ b/examples/AudioInpainting/Audio_Declipping_Example.m	Thu Jul 21 14:07:41 2011 +0100
@@ -33,32 +33,39 @@
 % subplot(1,2,2); imagesc(SMALL.Problem.imageSrc/SMALL.Problem.maxval);
 % title('Dictionary image:');colormap(gray);axis off; axis image;
 time=0;
+error2=0.001^2;
 coeffFrames = zeros(SMALL.Problem.p, SMALL.Problem.n);
 
 for i=1:SMALL.Problem.n
     
     idx = find(SMALL.Problem.M(:,i));
-    SMALL.Problem.A  = SMALL.Problem.B(idx,:);
+    if size(idx,1)==SMALL.Problem.m
+        continue
+    end
+    Dict = SMALL.Problem.B(idx,:);
+    wDict = 1./sqrt(diag(Dict'*Dict));
+    
+    SMALL.Problem.A  = Dict*diag(wDict);
     
     SMALL.Problem.b1 = SMALL.Problem.b(idx,i);
     
-   
+    
     
     %   Defining the parameters sparse representation
     SMALL.solver=SMALL_init_solver;
     SMALL.solver.toolbox='ompbox';
-    SMALL.solver.name='omp2';
+    SMALL.solver.name='omp2Gabor';
     
     SMALL.solver.param=struct(...
-        'epsilon', 0.001,...
-        'maxatoms', 64); 
+        'epsilon', error2*size(idx,1),...
+        'maxatoms', 128); 
     
     % Find solution
     
     SMALL.solver=SMALL_solve(SMALL.Problem, SMALL.solver);
     
     
-    coeffFrames(:,i) = SMALL.solver.solution;
+    coeffFrames(:,i) = wDict .* SMALL.solver.solution;
     time = time + SMALL.solver.time;
     
     
@@ -70,10 +77,32 @@
 SMALL.Problem.reconstruct=@(x) AudioDeclipping_reconstruct(x, SMALL.Problem);
 reconstructed=SMALL.Problem.reconstruct(coeffFrames);
 
-%%  plot time and psnr given dictionary size %%
-figure('Name', 'time and psnr');
 
-subplot(1,2,1); plot(dictsize(1,:), time(1,:), 'ro-');
-title('Time vs number of source image patches used');
-subplot(1,2,2); plot(dictsize(1,:), psnr(1,:), 'b*-');
-title('PSNR vs number of source image patches used');
\ No newline at end of file
+
+%% Plot results
+
+xClipped = SMALL.Problem.clipped;
+xClean   = SMALL.Problem.original;
+xEst1    = reconstructed.audioAllSamples;
+xEst2    = reconstructed.audioOnlyClipped;
+fs       = SMALL.Problem.fs;
+
+figure
+hold on
+plot(xClipped,'r')
+plot(xClean)
+plot(xEst2,'--g')
+plot([1;length(xClipped)],[1;1]*[-1,1]*max(abs(xClipped)),':r')
+legend('Clipped','True solution','Estimate')
+
+% Normalized and save sounds
+normX = 1.1*max(abs([xEst1(:);xEst2(:);xClean(:)]));
+L = min([length(xEst2),length(xEst1),length(xClean),length(xEst1),length(xClipped)]);
+xEst1 = xEst1(1:L)/normX;
+xEst2 = xEst2(1:L)/normX;
+xClipped = xClipped(1:L)/normX;
+xClean = xClean(1:L)/normX;
+wavwrite(xEst1,fs,[expParam.destDir 'xEst1.wav']);
+wavwrite(xEst2,fs,[expParam.destDir 'xEst2.wav']);
+wavwrite(xClipped,fs,[expParam.destDir 'xClipped.wav']);
+wavwrite(xClean,fs,[expParam.destDir 'xClean.wav']);
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/ksvd utils/ompbox utils/make.m	Thu Jul 21 14:07:41 2011 +0100
@@ -0,0 +1,41 @@
+function make
+%MAKE Build the OMPBox package.
+%  MAKE compiles all OMPBox MEX functions, using Matlab's default MEX
+%  compiler. If the MEX compiler has not been set-up before, please run
+%
+%    mex -setup
+%
+%  before using this MAKE file.
+
+%  Ron Rubinstein
+%  Computer Science Department
+%  Technion, Haifa 32000 Israel
+%  ronrubin@cs
+%
+%  August 2009
+
+
+% detect platform 
+
+compstr = computer;
+is64bit = strcmp(compstr(end-1:end),'64');
+
+
+% compilation parameters
+
+compile_params = cell(0);
+if (is64bit)
+  compile_params{1} = '-largeArrayDims';
+end
+
+
+% Compile files %
+
+ompsources = {'mexutils.c','ompcoreGabor.c','omputils.c','myblas.c','ompprof.c'};
+
+disp('Compiling ompmex...');
+mex('-g','ompmexGabor.c', ompsources{:},compile_params{:});
+
+disp('Compiling omp2mex...');
+mex('-g','omp2mexGabor.c',ompsources{:},compile_params{:});
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/ksvd utils/ompbox utils/mexutils.c	Thu Jul 21 14:07:41 2011 +0100
@@ -0,0 +1,79 @@
+/**************************************************************************
+ *
+ * File name: mexutils.c
+ *
+ * Ron Rubinstein
+ * Computer Science Department
+ * Technion, Haifa 32000 Israel
+ * ronrubin@cs
+ *
+ * Last Updated: 15.8.2009
+ *
+ *************************************************************************/
+
+#include "mexutils.h"
+#include <math.h>
+
+
+
+/* verify that the mxArray contains a double matrix */
+
+void checkmatrix(const mxArray *param, char *fname, char *pname)
+{
+  char errmsg[100];
+  sprintf(errmsg, "%.15s requires that %.25s be a double matrix.", fname, pname);
+  if (!mxIsDouble(param) || mxIsComplex(param) || mxGetNumberOfDimensions(param)>2) {
+    mexErrMsgTxt(errmsg);
+  }
+}
+
+
+/* verify that the mxArray contains a 1-D double vector */
+
+void checkvector(const mxArray *param, char *fname, char *pname)
+{
+  char errmsg[100];
+  sprintf(errmsg, "%.15s requires that %.25s be a double vector.", fname, pname);
+  if (!mxIsDouble(param) || mxIsComplex(param) || mxGetNumberOfDimensions(param)>2 || (mxGetM(param)!=1 && mxGetN(param)!=1)) {
+    mexErrMsgTxt(errmsg);
+  }
+}
+
+
+/* verify that the mxArray contains a double scalar */
+
+void checkscalar(const mxArray *param, char *fname, char *pname)
+{
+  char errmsg[100];
+  sprintf(errmsg, "%.15s requires that %.25s be a double scalar.", fname, pname);
+  if (!mxIsDouble(param) || mxIsComplex(param) || mxGetNumberOfDimensions(param)>2 || 
+      mxGetM(param)!=1 || mxGetN(param)!=1) 
+  {
+    mexErrMsgTxt(errmsg);
+  }
+}
+
+
+/* verify that the mxArray contains a sparse matrix */
+
+void checksparse(const mxArray *param, char *fname, char *pname)
+{
+  char errmsg[100];
+  sprintf(errmsg, "%.15s requires that %.25s be sparse.", fname, pname);
+  if (!mxIsSparse(param)) {
+    mexErrMsgTxt(errmsg);
+  }
+}
+
+
+/* verify that the mxArray contains a 1-dimensional cell array */
+
+void checkcell_1d(const mxArray *param, char *fname, char *pname)
+{
+  char errmsg[100];
+  sprintf(errmsg, "%.15s requires that %.25s be a 1-D cell array.", fname, pname);
+  if (!mxIsCell(param) || mxGetNumberOfDimensions(param)>2 || (mxGetM(param)!=1 && mxGetN(param)!=1)) {
+    mexErrMsgTxt(errmsg);
+  }
+}
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/ksvd utils/ompbox utils/mexutils.h	Thu Jul 21 14:07:41 2011 +0100
@@ -0,0 +1,103 @@
+/**************************************************************************
+ *
+ * File name: mexutils.h
+ *
+ * Ron Rubinstein
+ * Computer Science Department
+ * Technion, Haifa 32000 Israel
+ * ronrubin@cs
+ *
+ * Last Updated: 18.8.2009
+ *
+ * Utility functions for MEX files.
+ *
+ *************************************************************************/
+
+
+#ifndef __MEX_UTILS_H__
+#define __MEX_UTILS_H__
+
+#include "mex.h"
+
+
+
+/**************************************************************************
+ * Function checkmatrix:
+ *
+ * Verify that the specified mxArray is real, of type double, and has 
+ * no more than two dimensions. If not, an error message is printed
+ * and the mex file terminates.
+ * 
+ * Parameters:
+ *   param - the mxArray to be checked
+ *   fname - the name of the function where the error occured (15 characters or less)
+ *   pname - the name of the parameter (25 characters or less)
+ *
+ **************************************************************************/
+void checkmatrix(const mxArray *param, char *fname, char *pname);
+
+
+/**************************************************************************
+ * Function checkvector:
+ *
+ * Verify that the specified mxArray is 1-D, real, and of type double. The
+ * vector may be a column or row vector. Otherwise, an error message is
+ * printed and the mex file terminates.
+ * 
+ * Parameters:
+ *   param - the mxArray to be checked
+ *   fname - the name of the function where the error occured (15 characters or less)
+ *   pname - the name of the parameter (25 characters or less)
+ *
+ **************************************************************************/
+void checkvector(const mxArray *param, char *fname, char *pname);
+
+
+/**************************************************************************
+ * Function checkscalar:
+ *
+ * Verify that the specified mxArray represents a real double scalar value. 
+ * If not, an error message is printed and the mex file terminates.
+ * 
+ * Parameters:
+ *   param - the mxArray to be checked
+ *   fname - the name of the function where the error occured (15 characters or less)
+ *   pname - the name of the parameter (25 characters or less)
+ *
+ **************************************************************************/
+void checkscalar(const mxArray *param, char *fname, char *pname);
+
+
+/**************************************************************************
+ * Function checksparse:
+ *
+ * Verify that the specified mxArray contains a sparse matrix. If not,
+ * an error message is printed and the mex file terminates.
+ * 
+ * Parameters:
+ *   param - the mxArray to be checked
+ *   fname - the name of the function where the error occured (15 characters or less)
+ *   pname - the name of the parameter (25 characters or less)
+ *
+ **************************************************************************/
+void checksparse(const mxArray *param, char *fname, char *pname);
+
+
+/**************************************************************************
+ * Function checkcell_1d:
+ *
+ * Verify that the specified mxArray is a 1-D cell array. The cell array 
+ * may be arranged as either a column or a row. If not, an error message 
+ * is printed and the mex file terminates.
+ * 
+ * Parameters:
+ *   param - the mxArray to be checked
+ *   fname - the name of the function where the error occured (15 characters or less)
+ *   pname - the name of the parameter (25 characters or less)
+ *
+ **************************************************************************/
+void checkcell_1d(const mxArray *param, char *fname, char *pname);
+
+
+#endif
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/ksvd utils/ompbox utils/myblas.c	Thu Jul 21 14:07:41 2011 +0100
@@ -0,0 +1,673 @@
+/**************************************************************************
+ *
+ * File name: myblas.c
+ *
+ * Ron Rubinstein
+ * Computer Science Department
+ * Technion, Haifa 32000 Israel
+ * ronrubin@cs
+ *
+ * Version: 1.1
+ * Last updated: 13.8.2009
+ *
+ *************************************************************************/
+
+
+#include "myblas.h"
+#include <ctype.h>
+
+
+/* find maximum of absolute values */
+
+mwIndex maxabs(double c[], mwSize m)
+{
+  mwIndex maxid=0, k;
+  double absval, maxval = SQR(*c);   /* use square which is quicker than absolute value */
+
+  for (k=1; k<m; ++k) {
+    absval = SQR(c[k]);
+    if (absval > maxval) {
+      maxval = absval;
+      maxid = k;
+    }
+  }
+  return maxid;
+}
+
+
+/* compute y := alpha*x + y */
+
+void vec_sum(double alpha, double x[], double y[], mwSize n)
+{
+  mwIndex i;
+
+  for (i=0; i<n; ++i) {
+    y[i] += alpha*x[i];
+  }
+}
+
+/* compute y := alpha*x .* y */
+
+void vec_smult(double alpha, double x[], double y[], mwSize n)
+{
+  mwIndex i;
+
+  for (i=0; i<n; ++i) {
+    y[i] *= alpha*x[i];
+  }
+}
+
+/* compute y := alpha*A*x */
+
+void mat_vec(double alpha, double A[], double x[], double y[], mwSize n, mwSize m)
+{
+  mwIndex i, j, i_n;
+  double *Ax;
+
+  Ax = mxCalloc(n,sizeof(double));
+
+  for (i=0; i<m; ++i) {
+    i_n = i*n;
+    for (j=0; j<n; ++j) {
+      Ax[j] += A[i_n+j] * x[i];
+    }
+  }
+
+  for (j=0; j<n; ++j) {
+    y[j] = alpha*Ax[j];
+  }
+
+  mxFree(Ax);
+}
+
+
+/* compute y := alpha*A'*x */
+
+void matT_vec(double alpha, double A[], double x[], double y[], mwSize n, mwSize m)
+{
+  mwIndex i, j, n_i;
+  double sum0, sum1, sum2, sum3;
+
+  for (j=0; j<m; ++j) {
+    y[j] = 0;
+  }
+
+  /* use loop unrolling to accelerate computation */
+
+  for (i=0; i<m; ++i) {
+    n_i = n*i;
+    sum0 = sum1 = sum2 = sum3 = 0;
+    for (j=0; j+4<n; j+=4) {
+      sum0 += A[n_i+j]*x[j];
+      sum1 += A[n_i+j+1]*x[j+1];
+      sum2 += A[n_i+j+2]*x[j+2];
+      sum3 += A[n_i+j+3]*x[j+3];
+    }
+    y[i] += alpha * ((sum0 + sum1) + (sum2 + sum3));
+    while (j<n) {
+      y[i] += alpha*A[n_i+j]*x[j];
+      j++;
+    }
+  }
+}
+
+
+/* compute y := alpha*A*x */
+
+void mat_sp_vec(double alpha, double pr[], mwIndex ir[], mwIndex jc[], double x[], double y[], mwSize n, mwSize m)
+{
+  
+  mwIndex i, j, j1, j2;
+
+  for (i=0; i<n; ++i) {
+    y[i] = 0;
+  }
+  
+  j2 = jc[0];
+  for (i=0; i<m; ++i) {
+    j1 = j2; j2 = jc[i+1];
+    for (j=j1; j<j2; ++j) {
+      y[ir[j]] += alpha * pr[j] * x[i];
+    }
+  }
+  
+}
+
+
+/* compute y := alpha*A'*x */
+
+void matT_sp_vec(double alpha, double pr[], mwIndex ir[], mwIndex jc[], double x[], double y[], mwSize n, mwSize m)
+{
+  
+  mwIndex i, j, j1, j2;
+  
+  for (i=0; i<m; ++i) {
+    y[i] = 0;
+  }
+  
+  j2 = jc[0];
+  for (i=0; i<m; ++i) {
+    j1 = j2; j2 = jc[i+1];
+    for (j=j1; j<j2; ++j) {
+      y[i] += alpha * pr[j] * x[ir[j]];
+    }
+  }
+  
+}
+
+
+/* compute y := alpha*A*x */
+
+void mat_vec_sp(double alpha, double A[], double pr[], mwIndex ir[], mwIndex jc[], double y[], mwSize n, mwSize m)
+{
+  
+  mwIndex i, j, j_n, k, kend;
+  
+  for (i=0; i<n; ++i) {
+    y[i] = 0;
+  }
+  
+  kend = jc[1];
+  if (kend==0) {   /* x is empty */
+    return;
+  }
+  
+  for (k=0; k<kend; ++k) {
+    j = ir[k];
+    j_n = j*n;
+    for (i=0; i<n; ++i) {
+      y[i] += alpha * A[i+j_n] * pr[k];
+    }
+  }
+
+}
+
+
+/* compute y := alpha*A'*x */
+
+void matT_vec_sp(double alpha, double A[], double pr[], mwIndex ir[], mwIndex jc[], double y[], mwSize n, mwSize m)
+{
+  
+  mwIndex i, j, j_n, k, kend;
+  
+  for (i=0; i<m; ++i) {
+    y[i] = 0;
+  }
+  
+  kend = jc[1];
+  if (kend==0) {   /* x is empty */
+    return;
+  }
+  
+  for (j=0; j<m; ++j) {
+    j_n = j*n;
+    for (k=0; k<kend; ++k) {
+      i = ir[k];
+      y[j] += alpha * A[i+j_n] * pr[k];
+    }
+  }
+  
+}
+
+
+/* compute y := alpha*A*x */
+
+void mat_sp_vec_sp(double alpha, double pr[], mwIndex ir[], mwIndex jc[], double prx[], mwIndex irx[], mwIndex jcx[], double y[], mwSize n, mwSize m)
+{
+  
+  mwIndex i, j, k, kend, j1, j2;
+
+  for (i=0; i<n; ++i) {
+    y[i] = 0;
+  }
+  
+  kend = jcx[1]; 
+  if (kend==0) {   /* x is empty */
+    return;
+  }
+  
+  for (k=0; k<kend; ++k) {
+    i = irx[k];
+    j1 = jc[i]; j2 = jc[i+1];
+    for (j=j1; j<j2; ++j) {
+      y[ir[j]] += alpha * pr[j] * prx[k];
+    }
+  }
+  
+}
+
+
+/* compute y := alpha*A'*x */
+
+void matT_sp_vec_sp(double alpha, double pr[], mwIndex ir[], mwIndex jc[], double prx[], mwIndex irx[], mwIndex jcx[], double y[], mwSize n, mwSize m)
+{
+  
+  mwIndex i, j, k, jend, kend, jadd, kadd, delta;
+  
+  for (i=0; i<m; ++i) {
+    y[i] = 0;
+  }
+  
+  kend = jcx[1];
+  if (kend==0) {   /* x is empty */
+    return;
+  }
+  
+  for (i=0; i<m; ++i) {
+    j = jc[i]; 
+    jend = jc[i+1];
+    k = 0;
+    while (j<jend && k<kend) {
+      
+      delta = ir[j] - irx[k];
+      
+      if (delta) { /* if indices differ - increment the smaller one */
+        jadd = delta<0;
+        kadd = 1-jadd;
+        j += jadd;
+        k += kadd;
+      }
+      
+      else {    /* indices are equal - add to result and increment both */
+        y[i] += alpha * pr[j] * prx[k];
+        j++; k++;
+      }
+    }
+  }
+  
+}
+
+
+/* matrix-matrix multiplication */
+
+void mat_mat(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k)
+{
+  mwIndex i1, i2, i3, iX, iA, i2_n;
+  double b;
+  
+  for (i1=0; i1<n*k; i1++) {
+    X[i1] = 0;
+  }
+
+  for (i2=0; i2<m; ++i2) {
+    i2_n = i2*n;
+    iX = 0;
+    for (i3=0; i3<k; ++i3) {
+      iA = i2_n;
+      b = B[i2+i3*m];
+      for (i1=0; i1<n; ++i1) {
+        X[iX++] += A[iA++]*b;
+      }
+    }
+  }
+  
+  for (i1=0; i1<n*k; i1++) {
+    X[i1] *= alpha;
+  }
+}
+
+
+/* matrix-transpose-matrix multiplication */
+
+void matT_mat(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k)
+{
+  mwIndex i1, i2, i3, iX, iA, i2_n;
+  double *x, sum0, sum1, sum2, sum3;
+
+  for (i2=0; i2<m; ++i2) {
+    for (i3=0; i3<k; ++i3) {
+      sum0 = sum1 = sum2 = sum3 = 0;
+      for (i1=0; i1+4<n; i1+=4) {
+        sum0 += A[i1+0+i2*n]*B[i1+0+i3*n];
+        sum1 += A[i1+1+i2*n]*B[i1+1+i3*n];
+        sum2 += A[i1+2+i2*n]*B[i1+2+i3*n];
+        sum3 += A[i1+3+i2*n]*B[i1+3+i3*n];
+      }
+      X[i2+i3*m] = (sum0+sum1) + (sum2+sum3);
+      while(i1<n) {
+        X[i2+i3*m] += A[i1+i2*n]*B[i1+i3*n];
+        i1++;
+      }
+    }
+  }
+  
+  for (i1=0; i1<m*k; i1++) {
+    X[i1] *= alpha;
+  }
+}
+
+
+/* tensor-matrix product */
+
+void tens_mat(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k, mwSize l)
+{
+  mwIndex i1, i2, i3, i4, i2_n, nml;
+  double b;
+  
+  nml = n*m*l;
+  for (i1=0; i1<nml; ++i1) {
+    X[i1] = 0;
+  }
+
+  for (i2=0; i2<m; ++i2) {
+    i2_n = i2*n;
+    for (i3=0; i3<k; ++i3) {
+      for (i4=0; i4<l; ++i4) {
+        b = B[i4+i3*l];
+        for (i1=0; i1<n; ++i1) {
+          X[i1 + i2_n + i4*n*m] += A[i1 + i2_n + i3*n*m] * b;
+        }
+      }
+    }
+  }
+  
+  for (i1=0; i1<nml; ++i1) {
+    X[i1] *= alpha;
+  }
+}
+
+
+/* tensor-matrix-transpose product */
+
+void tens_matT(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k, mwSize l)
+{
+  mwIndex i1, i2, i3, i4, i2_n, nml;
+  double b;
+  
+  nml = n*m*l;
+  for (i1=0; i1<nml; ++i1) {
+    X[i1] = 0;
+  }
+
+  for (i2=0; i2<m; ++i2) {
+    i2_n = i2*n;
+    for (i4=0; i4<l; ++i4) {
+      for (i3=0; i3<k; ++i3) {
+        b = B[i3+i4*k];
+        for (i1=0; i1<n; ++i1) {
+          X[i1 + i2_n + i4*n*m] += A[i1 + i2_n + i3*n*m] * b;
+        }
+      }
+    }
+  }
+  
+  for (i1=0; i1<nml; ++i1) {
+    X[i1] *= alpha;
+  }
+}
+
+
+/* dot product */
+
+double dotprod(double a[], double b[], mwSize n)
+{
+  double sum = 0;
+  mwIndex i;
+  for (i=0; i<n; ++i)
+    sum += a[i]*b[i];
+  return sum;
+}
+
+
+/* find maximum of vector */
+
+mwIndex maxpos(double c[], mwSize m)
+{
+  mwIndex maxid=0, k;
+  double val, maxval = *c;
+
+  for (k=1; k<m; ++k) {
+    val = c[k];
+    if (val > maxval) {
+      maxval = val;
+      maxid = k;
+    }
+  }
+  return maxid;
+}
+
+
+/* solve L*x = b */
+
+void backsubst_L(double L[], double b[], double x[], mwSize n, mwSize k)
+{
+  mwIndex i, j;
+  double rhs;
+
+  for (i=0; i<k; ++i) {
+    rhs = b[i];
+    for (j=0; j<i; ++j) {
+      rhs -= L[j*n+i]*x[j];
+    }
+    x[i] = rhs/L[i*n+i];
+  }
+}
+
+
+/* solve L'*x = b */
+
+void backsubst_Lt(double L[], double b[], double x[], mwSize n, mwSize k)
+{
+  mwIndex i, j;
+  double rhs;
+
+  for (i=k; i>=1; --i) {
+    rhs = b[i-1];
+    for (j=i; j<k; ++j) {
+      rhs -= L[(i-1)*n+j]*x[j];
+    }
+    x[i-1] = rhs/L[(i-1)*n+i-1];
+  }
+}
+
+
+/* solve U*x = b */
+
+void backsubst_U(double U[], double b[], double x[], mwSize n, mwSize k)
+{
+  mwIndex i, j;
+  double rhs;
+
+  for (i=k; i>=1; --i) {
+    rhs = b[i-1];
+    for (j=i; j<k; ++j) {
+      rhs -= U[j*n+i-1]*x[j];
+    }
+    x[i-1] = rhs/U[(i-1)*n+i-1];
+  }
+}
+
+
+/* solve U'*x = b */
+
+void backsubst_Ut(double U[], double b[], double x[], mwSize n, mwSize k)
+{
+  mwIndex i, j;
+  double rhs;
+
+  for (i=0; i<k; ++i) {
+    rhs = b[i];
+    for (j=0; j<i; ++j) {
+      rhs -= U[i*n+j]*x[j];
+    }
+    x[i] = rhs/U[i*n+i];
+  }
+}
+
+
+/* back substitution solver */
+
+void backsubst(char ul, double A[], double b[], double x[], mwSize n, mwSize k)
+{
+  if (tolower(ul) == 'u') {
+    backsubst_U(A, b, x, n, k);
+  }
+  else if (tolower(ul) == 'l') {
+    backsubst_L(A, b, x, n, k);
+  }
+  else {
+    mexErrMsgTxt("Invalid triangular matrix type: must be ''U'' or ''L''");
+  }
+}
+
+
+/* solve equation set using cholesky decomposition */
+
+void cholsolve(char ul, double A[], double b[], double x[], mwSize n, mwSize k)
+{
+  double *tmp;
+
+  tmp = mxMalloc(k*sizeof(double));
+
+  if (tolower(ul) == 'l') {
+    backsubst_L(A, b, tmp, n, k);
+    backsubst_Lt(A, tmp, x, n, k);
+  }
+  else if (tolower(ul) == 'u') {
+    backsubst_Ut(A, b, tmp, n, k);
+    backsubst_U(A, tmp, x, n, k);
+  }
+  else {
+    mexErrMsgTxt("Invalid triangular matrix type: must be either ''U'' or ''L''");
+  }
+
+  mxFree(tmp);
+}
+
+
+/* perform a permutation assignment y := x(ind(1:k)) */
+
+void vec_assign(double y[], double x[], mwIndex ind[], mwSize k)
+{
+  mwIndex i;
+
+  for (i=0; i<k; ++i)
+    y[i] = x[ind[i]];
+}
+
+
+/* matrix transpose */
+
+void transpose(double X[], double Y[], mwSize n, mwSize m)
+{
+  mwIndex i, j, i_m, j_n;
+  
+  if (n<m) {
+    for (j=0; j<m; ++j) {
+      j_n = j*n;
+      for (i=0; i<n; ++i) {
+        Y[j+i*m] = X[i+j_n];
+      }
+    }
+  }
+  else {
+    for (i=0; i<n; ++i) {
+      i_m = i*m;
+      for (j=0; j<m; ++j) {
+        Y[j+i_m] = X[i+j*n];
+      }
+    }
+  }
+}
+
+
+/* print contents of matrix */
+
+void printmat(double A[], int n, int m, char* matname)
+{
+  int i, j;
+  mexPrintf("\n%s = \n\n", matname);
+
+  if (n*m==0) {
+    mexPrintf("   Empty matrix: %d-by-%d\n\n", n, m);
+    return;
+  }
+
+  for (i=0; i<n; ++i) {
+    for (j=0; j<m; ++j)
+      mexPrintf("   %lf", A[j*n+i]);
+    mexPrintf("\n");
+  }
+  mexPrintf("\n");
+}
+
+
+/* print contents of sparse matrix */
+
+void printspmat(mxArray *a, char* matname)
+{
+  mwIndex *aJc = mxGetJc(a);
+  mwIndex *aIr = mxGetIr(a);
+  double *aPr = mxGetPr(a);
+
+  int i;
+
+  mexPrintf("\n%s = \n\n", matname);
+
+  for (i=0; i<aJc[1]; ++i)
+    printf("   (%d,1) = %lf\n", aIr[i]+1,aPr[i]);
+
+  mexPrintf("\n");
+}
+
+
+
+/* matrix multiplication using Winograd's algorithm */
+
+/*
+void mat_mat2(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k)
+{
+  
+  mwIndex i1, i2, i3, iX, iA, i2_n;
+  double b, *AA, *BB;
+  
+  AA = mxCalloc(n,sizeof(double));
+  BB = mxCalloc(k,sizeof(double));
+  
+  for (i1=0; i1<n*k; i1++) {
+    X[i1] = 0;
+  }
+  
+  for (i1=0; i1<n; ++i1) {
+    for (i2=0; i2<m/2; ++i2) {
+      AA[i1] += A[i1+2*i2*n]*A[i1+(2*i2+1)*n];
+    }
+  }
+
+  for (i2=0; i2<k; ++i2) {
+    for (i1=0; i1<m/2; ++i1) {
+      BB[i2] += B[2*i1+i2*m]*B[2*i1+1+i2*m];
+    }
+  }
+
+  for (i2=0; i2<k; ++i2) {
+    for (i3=0; i3<m/2; ++i3) {
+      for (i1=0; i1<n; ++i1) {
+        X[i1+i2*n] += (A[i1+(2*i3)*n]+B[2*i3+1+i2*m])*(A[i1+(2*i3+1)*n]+B[2*i3+i2*m]);
+      }
+    }
+  }
+  
+  if (m%2) {
+    for (i2=0; i2<k; ++i2) {
+      for (i1=0; i1<n; ++i1) {
+        X[i1+i2*n] += A[i1+(m-1)*n]*B[m-1+i2*m];
+      }
+    }
+  }
+  
+  for (i2=0; i2<k; ++i2) {
+    for (i1=0; i1<n; ++i1) {
+      X[i1+i2*n] -= (AA[i1] + BB[i2]);
+      X[i1+i2*n] *= alpha;
+    }
+  }
+  
+  mxFree(AA);
+  mxFree(BB);
+}
+*/
+
+
+
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/ksvd utils/ompbox utils/myblas.h	Thu Jul 21 14:07:41 2011 +0100
@@ -0,0 +1,511 @@
+/**************************************************************************
+ *
+ * File name: myblas.h
+ *
+ * Ron Rubinstein
+ * Computer Science Department
+ * Technion, Haifa 32000 Israel
+ * ronrubin@cs
+ *
+ * Version: 1.1
+ * Last updated: 17.8.2009
+ *
+ * A collection of basic linear algebra functions, in the spirit of the
+ * BLAS/LAPACK libraries.
+ *
+ *************************************************************************/
+
+
+
+#ifndef __MY_BLAS_H__
+#define __MY_BLAS_H__
+
+
+#include "mex.h"
+#include <math.h>
+
+
+
+/**************************************************************************
+ * Squared value.
+ **************************************************************************/
+#define SQR(X) ((X)*(X))
+
+
+
+/**************************************************************************
+ * Matrix-vector multiplication. 
+ *
+ * Computes an operation of the form:
+ *
+ *   y := alpha*A*x
+ *
+ * Parameters:
+ *   A - matrix of size n X m
+ *   x - vector of length m
+ *   y - output vector of length n
+ *   alpha - real constant
+ *   n, m - dimensions of A
+ *
+ * Note: This function re-writes the contents of y.
+ *
+ **************************************************************************/
+void mat_vec(double alpha, double A[], double x[], double y[], mwSize n, mwSize m);
+
+
+
+/**************************************************************************
+ * Matrix-transpose-vector multiplication. 
+ *
+ * Computes an operation of the form:
+ *
+ *   y := alpha*A'*x
+ *
+ * Parameters:
+ *   A - matrix of size n X m
+ *   x - vector of length n
+ *   y - output vector of length m
+ *   alpha - real constant
+ *   n, m - dimensions of A
+ *
+ * Note: This function re-writes the contents of y.
+ *
+ **************************************************************************/
+void matT_vec(double alpha, double A[], double x[], double y[], mwSize n, mwSize m);
+
+
+
+/**************************************************************************
+ * Sparse-matrix-vector multiplication. 
+ *
+ * Computes an operation of the form:
+ *
+ *   y := alpha*A*x
+ *
+ * where A is a sparse matrix.
+ *
+ * Parameters:
+ *   pr,ir,jc - sparse representation of the matrix A, of size n x m
+ *   x - vector of length m
+ *   y - output vector of length n
+ *   alpha - real constant
+ *   n, m - dimensions of A
+ *
+ * Note: This function re-writes the contents of y.
+ *
+ **************************************************************************/
+void mat_sp_vec(double alpha, double pr[], mwIndex ir[], mwIndex jc[], double x[], double y[], mwSize n, mwSize m);
+
+
+
+/**************************************************************************
+ * Sparse-matrix-transpose-vector multiplication. 
+ *
+ * Computes an operation of the form:
+ *
+ *   y := alpha*A'*x
+ *
+ * where A is a sparse matrix.
+ *
+ * Parameters:
+ *   pr,ir,jc - sparse representation of the matrix A, of size n x m
+ *   x - vector of length m
+ *   y - output vector of length n
+ *   alpha - real constant
+ *   n, m - dimensions of A
+ *
+ * Note: This function re-writes the contents of y.
+ *
+ **************************************************************************/
+void matT_sp_vec(double alpha, double pr[], mwIndex ir[], mwIndex jc[], double x[], double y[], mwSize n, mwSize m);
+
+
+
+/**************************************************************************
+ * Matrix-sparse-vector multiplication. 
+ *
+ * Computes an operation of the form:
+ *
+ *   y := alpha*A*x
+ *
+ * where A is a matrix and x is a sparse vector.
+ *
+ * Parameters:
+ *   A - matrix of size n X m
+ *   pr,ir,jc - sparse representation of the vector x, of length m
+ *   y - output vector of length n
+ *   alpha - real constant
+ *   n, m - dimensions of A
+ *
+ * Note: This function re-writes the contents of y.
+ *
+ **************************************************************************/
+void mat_vec_sp(double alpha, double A[], double pr[], mwIndex ir[], mwIndex jc[], double y[], mwSize n, mwSize m);
+
+
+
+/**************************************************************************
+ * Matrix-transpose-sparse-vector multiplication. 
+ *
+ * Computes an operation of the form:
+ *
+ *   y := alpha*A'*x
+ *
+ * where A is a matrix and x is a sparse vector.
+ *
+ * Parameters:
+ *   A - matrix of size n X m
+ *   pr,ir,jc - sparse representation of the vector x, of length n
+ *   y - output vector of length m
+ *   alpha - real constant
+ *   n, m - dimensions of A
+ *
+ * Note: This function re-writes the contents of y.
+ *
+ **************************************************************************/
+void matT_vec_sp(double alpha, double A[], double pr[], mwIndex ir[], mwIndex jc[], double y[], mwSize n, mwSize m);
+
+
+
+/**************************************************************************
+ * Sparse-matrix-sparse-vector multiplication. 
+ *
+ * Computes an operation of the form:
+ *
+ *   y := alpha*A*x
+ *
+ * where A is a sparse matrix and x is a sparse vector.
+ *
+ * Parameters:
+ *   pr,ir,jc - sparse representation of the matrix A, of size n x m
+ *   prx,irx,jcx - sparse representation of the vector x (of length m)
+ *   y - output vector of length n
+ *   alpha - real constant
+ *   n, m - dimensions of A
+ *
+ * Note: This function re-writes the contents of y.
+ *
+ **************************************************************************/
+void mat_sp_vec_sp(double alpha, double pr[], mwIndex ir[], mwIndex jc[], double prx[], mwIndex irx[], mwIndex jcx[], double y[], mwSize n, mwSize m);
+
+
+
+/**************************************************************************
+ * Sparse-matrix-transpose-sparse-vector multiplication. 
+ *
+ * Computes an operation of the form:
+ *
+ *   y := alpha*A'*x
+ *
+ * where A is a sparse matrix and x is a sparse vector.
+ *
+ * Importnant note: this function is provided for completeness, but is NOT efficient.
+ * If possible, convert x to non-sparse representation and use matT_vec_sp instead.
+ *
+ * Parameters:
+ *   pr,ir,jc - sparse representation of the matrix A, of size n x m
+ *   prx,irx,jcx - sparse representation of the vector x (of length n)
+ *   y - output vector of length n
+ *   alpha - real constant
+ *   n, m - dimensions of A
+ *
+ * Note: This function re-writes the contents of y.
+ *
+ **************************************************************************/
+void matT_sp_vec_sp(double alpha, double pr[], mwIndex ir[], mwIndex jc[], double prx[], mwIndex irx[], mwIndex jcx[], double y[], mwSize n, mwSize m);
+
+
+
+/**************************************************************************
+ * Matrix-matrix multiplication. 
+ *
+ * Computes an operation of the form:
+ *
+ *   X := alpha*A*B
+ *
+ * Parameters:
+ *   A - matrix of size n X m
+ *   B - matrix of size m X k
+ *   X - output matrix of size n X k
+ *   alpha - real constant
+ *   n, m, k - dimensions of A, B
+ *
+ * Note: This function re-writes the contents of X.
+ *
+ **************************************************************************/
+void mat_mat(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k);
+
+
+
+/**************************************************************************
+ * Matrix-transpose-matrix multiplication. 
+ *
+ * Computes an operation of the form:
+ *
+ *   X := alpha*A*B
+ *
+ * Parameters:
+ *   A - matrix of size n X m
+ *   B - matrix of size m X k
+ *   X - output matrix of size n X k
+ *   alpha - real constant
+ *   n, m, k - dimensions of A, B
+ *
+ * Note: This function re-writes the contents of X.
+ *
+ **************************************************************************/
+void matT_mat(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k);
+
+
+
+/**************************************************************************
+ * Tensor-matrix multiplication. 
+ *
+ * This function accepts a 3-D tensor A of size n X m X k
+ * and a 2-D matrix B of size l X k.
+ * The function computes the 3-D tensor X of size n X m X l, where
+ *
+ *   X(i,j,:) = B*A(i,j,:)
+ *
+ * for all i,j.
+ *
+ * Parameters:
+ *   A - tensor of size n X m X k
+ *   B - matrix of size l X k
+ *   X - output tensor of size n X m X l
+ *   alpha - real constant
+ *   n, m, k, l - dimensions of A, B
+ *
+ * Note: This function re-writes the contents of X.
+ *
+ **************************************************************************/
+void tens_mat(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k, mwSize l);
+
+
+
+/**************************************************************************
+ * Tensor-matrix-transpose multiplication. 
+ *
+ * This function accepts a 3-D tensor A of size n X m X k
+ * and a 2-D matrix B of size k X l.
+ * The function computes the 3-D tensor X of size n X m X l, where
+ *
+ *   X(i,j,:) = B'*A(i,j,:)
+ *
+ * for all i,j.
+ *
+ * Parameters:
+ *   A - tensor of size n X m X k
+ *   B - matrix of size k X l
+ *   X - output tensor of size n X m X l
+ *   alpha - real constant
+ *   n, m, k, l - dimensions of A, B
+ *
+ * Note: This function re-writes the contents of X.
+ *
+ **************************************************************************/
+void tens_matT(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k, mwSize l);
+
+
+
+/**************************************************************************
+ * Vector-vector sum.
+ *
+ * Computes an operation of the form:
+ *
+ *   y := alpha*x + y
+ *
+ * Parameters:
+ *   x - vector of length n
+ *   y - output vector of length n
+ *   alpha - real constant
+ *   n - length of x,y
+ *
+ * Note: This function re-writes the contents of y.
+ *
+ **************************************************************************/
+void vec_sum(double alpha, double x[], double y[], mwSize n);
+
+/**************************************************************************
+ * Vector-vector scalar multiply.
+ *
+ * Computes an operation of the form:
+ *
+ *   y := alpha* x.*y
+ *
+ * Parameters:
+ *   x - vector of length n
+ *   y - output vector of length n
+ *   alpha - real constant
+ *   n - length of x,y
+ *
+ * Note: This function re-writes the contents of y.
+ *
+ **************************************************************************/
+
+
+void vec_smult(double alpha, double x[], double y[], mwSize n);
+
+
+/**************************************************************************
+ * Triangular back substitution.
+ *
+ * Solve the set of linear equations
+ *
+ *   T*x = b
+ *
+ * where T is lower or upper triangular.
+ *
+ * Parameters:
+ *   ul - 'U' for upper triangular, 'L' for lower triangular
+ *   A  - matrix of size n x m containing T
+ *   b  - vector of length k
+ *   x  - output vector of length k
+ *   n  - size of first dimension of A
+ *   k  - the size of the equation set, k<=n,m
+ *
+ * Note:
+ *   The matrix A can be of any size n X m, as long as n,m >= k. 
+ *   Only the lower/upper triangle of the submatrix A(1:k,1:k) defines the
+ *   matrix T (depending on the parameter ul).
+ *
+ **************************************************************************/
+void backsubst(char ul, double A[], double b[], double x[], mwSize n, mwSize k);
+
+
+
+/**************************************************************************
+ * Solve a set of equations using a Cholesky decomposition.
+ *
+ * Solve the set of linear equations
+ *
+ *   M*x = b
+ *
+ * where M is positive definite with a known Cholesky decomposition:
+ * either M=L*L' (L lower triangular) or M=U'*U (U upper triangular).
+ *
+ * Parameters:
+ *   ul - 'U' for upper triangular, 'L' for lower triangular decomposition
+ *   A  - matrix of size n x m with the Cholesky decomposition of M
+ *   b  - vector of length k
+ *   x  - output vector of length k
+ *   n  - size of first dimension of A
+ *   k  - the size of the equation set, k<=n,m
+ *
+ * Note:
+ *   The matrix A can be of any size n X m, as long as n,m >= k. 
+ *   Only the lower/upper triangle of the submatrix A(1:k,1:k) is used as
+ *   the Cholesky decomposition of M (depending on the parameter ul).
+ *
+ **************************************************************************/
+void cholsolve(char ul, double A[], double b[], double x[], mwSize n, mwSize k);
+
+
+
+/**************************************************************************
+ * Maximum absolute value.
+ *
+ * Returns the index of the coefficient with maximal absolute value in a vector.
+ *
+ * Parameters:
+ *   x - vector of length n
+ *   n - length of x
+ *
+ **************************************************************************/
+mwIndex maxabs(double x[], mwSize n);
+
+
+
+/**************************************************************************
+ * Maximum vector element.
+ *
+ * Returns the index of the maximal coefficient in a vector.
+ *
+ * Parameters:
+ *   x - vector of length n
+ *   n - length of x
+ *
+ **************************************************************************/
+mwIndex maxpos(double x[], mwSize n);
+
+
+
+/**************************************************************************
+ * Vector-vector dot product.
+ *
+ * Computes an operation of the form:
+ *
+ *   c = a'*b
+ *
+ * Parameters:
+ *   a, b - vectors of length n
+ *   n - length of a,b
+ *
+ * Returns: The dot product c.
+ *
+ **************************************************************************/
+double dotprod(double a[], double b[], mwSize n);
+
+
+
+/**************************************************************************
+ * Indexed vector assignment.
+ *
+ * Perform a permutation assignment of the form
+ *
+ *   y = x(ind)
+ *
+ * where ind is an array of indices to x.
+ *
+ * Parameters:
+ *   y - output vector of length k
+ *   x - input vector of arbitrary length
+ *   ind - array of indices into x (indices begin at 0)
+ *   k - length of the array ind
+ *
+ **************************************************************************/
+void vec_assign(double y[], double x[], mwIndex ind[], mwSize k);
+
+
+
+/**************************************************************************
+ * Matrix transpose.
+ *
+ * Computes Y := X'
+ *
+ * Parameters:
+ *   X - input matrix of size n X m
+ *   Y - output matrix of size m X n
+ *   n, m - dimensions of X
+ *
+ **************************************************************************/
+void transpose(double X[], double Y[], mwSize n, mwSize m);
+
+
+
+/**************************************************************************
+ * Print a matrix.
+ *
+ * Parameters:
+ *   A - matrix of size n X m
+ *   n, m - dimensions of A
+ *   matname - name of matrix to display
+ *
+ **************************************************************************/
+void printmat(double A[], int n, int m, char* matname);
+
+
+
+/**************************************************************************
+ * Print a sparse matrix.
+ *
+ * Parameters:
+ *   A - sparse matrix of type double
+ *   matname - name of matrix to display
+ *
+ **************************************************************************/
+void printspmat(mxArray *A, char* matname);
+
+
+#endif
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/ksvd utils/ompbox utils/omp2Gabor.m	Thu Jul 21 14:07:41 2011 +0100
@@ -0,0 +1,205 @@
+function gamma = omp2(varargin)
+%OMP2 Error-constrained Orthogonal Matching Pursuit.
+%  GAMMA = OMP2(D,X,G,EPSILON) solves the optimization problem
+%
+%       min  |GAMMA|_0     s.t.  |X - D*GAMMA|_2 <= EPSILON
+%      gamma
+%
+%  for each of the signals in X, using Batch Orthogonal Matching Pursuit.
+%  Here, D is a dictionary with normalized columns, X is a matrix
+%  containing column signals, EPSILON is the error target for each signal,
+%  and G is the Gramm matrix D'*D. The output GAMMA is a matrix containing
+%  the sparse representations as its columns. 
+%
+%  GAMMA = OMP2(D,X,[],EPSILON) performs the same operation, but without
+%  the matrix G, using OMP-Cholesky. This call produces the same output as
+%  Batch-OMP, but is significantly slower. Using this syntax is only
+%  recommended when available memory is too small to store G.
+%
+%  GAMMA = OMP2(DtX,XtX,G,EPSILON) is the fastest implementation of OMP2,
+%  but also requires the most memory. Here, DtX stores the projections
+%  D'*X, and XtX is a row vector containing the squared norms of the
+%  signals, sum(X.*X). In this case Batch-OMP is used, but without having
+%  to compute D'*X and XtX in advance, which slightly improves runtime.
+%  Note that in general, the call
+%
+%    GAMMA = OMP2(D'*X, sum(X.*X), G, EPSILON);
+%
+%  will be faster than the call
+%
+%    GAMMA = OMP2(D,X,G,EPSILON);
+%
+%  due to optimized matrix multiplications in Matlab. However, when the
+%  entire matrix D'*X cannot be stored in memory, one of the other two
+%  versions can be used. Both compute D'*X for just one signal at a time,
+%  and thus require much less memory.
+%
+%  GAMMA = OMP2(...,PARAM1,VAL1,PARAM2,VAL2,...) specifies additional
+%  parameters for OMP2. Available parameters are:
+%
+%    'gammamode' - Specifies the representation mode for GAMMA. Can be
+%                  either 'full' or 'sparse', corresponding to a full or
+%                  sparse matrix, respectively. By default, GAMMA is
+%                  returned as a sparse matrix.
+%    'maxatoms' -  Limits the number of atoms in the representation of each
+%                  signal. If specified, the number of atoms in each
+%                  representation does not exceed this number, even if the
+%                  error target is not met. Specifying maxatoms<0 implies
+%                  no limit (default).
+%    'messages'  - Specifies whether progress messages should be displayed.
+%                  When positive, this is the number of seconds between
+%                  status prints. When negative, indicates that no messages
+%                  should be displayed (this is the default).
+%    'checkdict' - Specifies whether dictionary normalization should be
+%                  verified. When set to 'on' (default) the dictionary
+%                  atoms are verified to be of unit L2-norm. Setting this
+%                  parameter to 'off' disables verification and accelerates
+%                  function performance. Note that an unnormalized
+%                  dictionary will produce invalid results.
+%    'profile'   - Can be either 'on' or 'off'. When 'on', profiling
+%                  information is displayed at the end of the funciton
+%                  execution.
+%
+%
+%  Summary of OMP2 versions:
+%
+%    version                 |   speed     |   memory
+%  -------------------------------------------------------------
+%   OMP2(DtX,XtX,G,EPSILON)  |  very fast  |  very large
+%   OMP2(D,X,G,EPSILON)      |  fast       |  moderate
+%   OMP2(D,X,[],EPSILON)     |  very slow  |  small
+%  -------------------------------------------------------------
+%
+%
+%  References:
+%  [1] M. Elad, R. Rubinstein, and M. Zibulevsky, "Efficient Implementation
+%      of the K-SVD Algorithm using Batch Orthogonal Matching Pursuit",
+%      Technical Report - CS, Technion, April 2008.
+%
+%  See also OMP.
+
+
+%  Ron Rubinstein
+%  Computer Science Department
+%  Technion, Haifa 32000 Israel
+%  ronrubin@cs
+%
+%  April 2009
+
+
+% default options
+
+sparse_gamma = 1;
+msgdelta = -1;
+maxatoms = -1;
+checkdict = 1;
+profile = 0;
+
+
+% determine number of parameters
+
+paramnum = 1;
+while (paramnum<=nargin && ~ischar(varargin{paramnum}))
+  paramnum = paramnum+1;
+end
+paramnum = paramnum-1;
+
+
+% parse options
+
+for i = paramnum+1:2:length(varargin)
+  paramname = varargin{i};
+  paramval = varargin{i+1};
+
+  switch lower(paramname)
+
+    case 'gammamode'
+      if (strcmpi(paramval,'sparse'))
+        sparse_gamma = 1;
+      elseif (strcmpi(paramval,'full'))
+        sparse_gamma = 0;
+      else
+        error('Invalid GAMMA mode');
+      end
+      
+    case 'maxatoms'
+      maxatoms = paramval;
+      
+    case 'messages'
+      msgdelta = paramval;
+
+    case 'checkdict'
+      if (strcmpi(paramval,'on'))
+        checkdict = 1;
+      elseif (strcmpi(paramval,'off'))
+        checkdict = 0;
+      else
+        error('Invalid checkdict option');
+      end
+
+    case 'profile'
+      if (strcmpi(paramval,'on'))
+        profile = 1;
+      elseif (strcmpi(paramval,'off'))
+        profile = 0;
+      else
+        error('Invalid profile mode');
+      end
+
+    otherwise
+      error(['Unknown option: ' paramname]);
+  end
+  
+end
+
+
+% determine call type
+
+if (paramnum==4)
+  
+  n1 = size(varargin{1},1);
+  n2 = size(varargin{2},1);
+  n3 = size(varargin{3},1);
+  
+  if ( (n1>1 && n2==1) || (n1==1 && n2==1 && n3==1) )  %  DtX,XtX,G,EPSILON
+    
+    DtX = varargin{1};
+    XtX = varargin{2};
+    G = varargin{3};
+    epsilon = varargin{4};
+    D = [];
+    X = [];
+    
+  else  % D,X,G,EPSILON
+    
+    D = varargin{1};
+    X = varargin{2};
+    G = varargin{3};
+    epsilon = varargin{4};
+    DtX = [];
+    XtX = [];
+    
+  end
+  
+else
+  error('Invalid number of parameters');
+end
+
+
+% verify dictionary normalization
+
+if (checkdict)
+  if (isempty(G))
+    atomnorms = sum(D.*D);
+  else
+    atomnorms = diag(G);
+  end
+  if (any(abs(atomnorms-1) > 1e-2))
+    error('Dictionary columns must be normalized to unit length');
+  end
+end
+
+
+% omp
+
+gamma = omp2mexGabor(D,X,DtX,XtX,G,epsilon,sparse_gamma,msgdelta,maxatoms,profile);
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/ksvd utils/ompbox utils/omp2mex.c	Thu Jul 21 14:07:41 2011 +0100
@@ -0,0 +1,156 @@
+/**************************************************************************
+ *
+ * File name: omp2mex.c
+ *
+ * Ron Rubinstein
+ * Computer Science Department
+ * Technion, Haifa 32000 Israel
+ * ronrubin@cs
+ *
+ * Last Updated: 18.8.2009
+ *
+ *************************************************************************/
+
+#include "ompcore.h"
+#include "omputils.h"
+#include "mexutils.h"
+
+
+/* Input Arguments */
+
+#define	IN_D	        prhs[0]
+#define IN_X          prhs[1]
+#define IN_DtX        prhs[2]
+#define IN_XtX        prhs[3]
+#define IN_G          prhs[4]
+#define IN_EPS        prhs[5]
+#define IN_SPARSE_G   prhs[6]
+#define IN_MSGDELTA   prhs[7]
+#define IN_MAXATOMS   prhs[8]
+#define IN_PROFILE    prhs[9]
+
+
+/* Output Arguments */
+
+#define	GAMMA_OUT     plhs[0]
+
+
+/***************************************************************************************/
+
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray*prhs[])
+
+{
+  double *D, *x, *DtX, *XtX, *G, eps, msgdelta;
+  int gmode, maxatoms, profile;
+  mwSize m, n, L;    /* D is n x m , X is n x L, DtX is m x L */
+
+  
+  /* check parameters */
+  
+  checkmatrix(IN_D, "OMP2", "D");
+  checkmatrix(IN_X, "OMP2", "X");
+  checkmatrix(IN_DtX, "OMP2", "DtX");
+  checkmatrix(IN_XtX, "OMP2", "XtX");
+  checkmatrix(IN_G, "OMP2", "G");
+  
+  checkscalar(IN_EPS, "OMP2", "EPSILON");
+  checkscalar(IN_SPARSE_G, "OMP2", "sparse_g");
+  checkscalar(IN_MSGDELTA, "OMP2", "msgdelta");
+  checkscalar(IN_MAXATOMS, "OMP2", "maxatoms");
+  checkscalar(IN_PROFILE, "OMP2", "profile");
+  
+  
+  /* get parameters */
+  
+  x = D = DtX = XtX = G = 0;
+  
+  if (!mxIsEmpty(IN_D))
+    D = mxGetPr(IN_D);
+  
+  if (!mxIsEmpty(IN_X))
+    x = mxGetPr(IN_X);
+  
+  if (!mxIsEmpty(IN_DtX))
+    DtX = mxGetPr(IN_DtX);
+  
+  if (!mxIsEmpty(IN_XtX))
+    XtX = mxGetPr(IN_XtX);
+  
+  if (!mxIsEmpty(IN_G))
+    G = mxGetPr(IN_G);
+  
+  eps = mxGetScalar(IN_EPS);
+  if ((int)(mxGetScalar(IN_SPARSE_G)+1e-2)) {
+    gmode = SPARSE_GAMMA;
+  }
+  else {
+    gmode = FULL_GAMMA;
+  }
+  msgdelta = mxGetScalar(IN_MSGDELTA);
+  if (mxGetScalar(IN_MAXATOMS) < -1e-5) {
+    maxatoms = -1;
+  }
+  else {
+    maxatoms = (int)(mxGetScalar(IN_MAXATOMS)+1e-2);
+  }
+  profile = (int)(mxGetScalar(IN_PROFILE)+1e-2);
+  
+  
+  /* check sizes */
+  
+  if (D && x) {
+    n = mxGetM(IN_D);
+    m = mxGetN(IN_D);
+    L = mxGetN(IN_X);
+    
+    if (mxGetM(IN_X) != n) {
+      mexErrMsgTxt("D and X have incompatible sizes.");
+    }
+    
+    if (G) {
+      if (mxGetN(IN_G)!=mxGetM(IN_G)) {
+        mexErrMsgTxt("G must be a square matrix.");
+      }
+      if (mxGetN(IN_G) != m) {
+        mexErrMsgTxt("D and G have incompatible sizes.");
+      }
+    }
+  }
+  
+  else if (DtX && XtX) {
+    m = mxGetM(IN_DtX);
+    L = mxGetN(IN_DtX);
+    
+    /* set n to an arbitrary value that is at least the max possible number of selected atoms */
+    
+    if (maxatoms>0) {
+      n = maxatoms;
+    }
+    else {
+      n = m;
+    }
+    
+    if ( !(mxGetM(IN_XtX)==L && mxGetN(IN_XtX)==1) && !(mxGetM(IN_XtX)==1 && mxGetN(IN_XtX)==L) ) {
+      mexErrMsgTxt("DtX and XtX have incompatible sizes.");
+    }
+    
+    if (mxGetN(IN_G)!=mxGetM(IN_G)) {
+      mexErrMsgTxt("G must be a square matrix.");
+    }
+    if (mxGetN(IN_G) != m) {
+      mexErrMsgTxt("DtX and G have incompatible sizes.");
+    }
+  }
+  
+  else {
+    mexErrMsgTxt("Either D and X, or DtX and XtX, must be specified.");
+  }
+  
+  
+  /* Do OMP! */
+  
+  GAMMA_OUT = ompcore(D, x, DtX, XtX, G, n, m, L, maxatoms, eps, gmode, profile, msgdelta, 1);
+  
+  return;
+}
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/ksvd utils/ompbox utils/omp2mex.m	Thu Jul 21 14:07:41 2011 +0100
@@ -0,0 +1,23 @@
+%This is the Matlab interface to the OMP2 MEX implementation.
+%The function is not for independent use, only through omp2.m.
+
+
+%OMP2MEX Matlab interface to the OMP2 MEX implementation.
+%  GAMMA = OMP2MEX(D,X,DtX,XtX,G,EPSILON,SPARSE_G,MSGDELTA,MAXATOMS,PROFILE)
+%  invokes the OMP2 MEX function according to the specified parameters. Not
+%  all the parameters are required. Those among D, X, DtX, XtX and G which
+%  are not specified should be passed as [].
+%
+%  EPSILON - the target error.
+%  SPARSE_G - returns a sparse GAMMA when nonzero, full GAMMA when zero.
+%  MSGDELTA - the delay in secs between messages. Zero means no messages.
+%  MAXATOMS - the max number of atoms per signal, negative for no max.
+%  PROFILE - nonzero means that profiling information should be printed.
+
+
+%  Ron Rubinstein
+%  Computer Science Department
+%  Technion, Haifa 32000 Israel
+%  ronrubin@cs
+%
+%  April 2009
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/ksvd utils/ompbox utils/omp2mexGabor.c	Thu Jul 21 14:07:41 2011 +0100
@@ -0,0 +1,156 @@
+/**************************************************************************
+ *
+ * File name: omp2mex.c
+ *
+ * Ron Rubinstein
+ * Computer Science Department
+ * Technion, Haifa 32000 Israel
+ * ronrubin@cs
+ *
+ * Last Updated: 18.8.2009
+ *
+ *************************************************************************/
+
+#include "ompcoreGabor.h"
+#include "omputils.h"
+#include "mexutils.h"
+
+
+/* Input Arguments */
+
+#define	IN_D	        prhs[0]
+#define IN_X          prhs[1]
+#define IN_DtX        prhs[2]
+#define IN_XtX        prhs[3]
+#define IN_G          prhs[4]
+#define IN_EPS        prhs[5]
+#define IN_SPARSE_G   prhs[6]
+#define IN_MSGDELTA   prhs[7]
+#define IN_MAXATOMS   prhs[8]
+#define IN_PROFILE    prhs[9]
+
+
+/* Output Arguments */
+
+#define	GAMMA_OUT     plhs[0]
+
+
+/***************************************************************************************/
+
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray*prhs[])
+
+{
+  double *D, *x, *DtX, *XtX, *G, eps, msgdelta;
+  int gmode, maxatoms, profile;
+  mwSize m, n, L;    /* D is n x m , X is n x L, DtX is m x L */
+
+  
+  /* check parameters */
+  
+  checkmatrix(IN_D, "OMP2", "D");
+  checkmatrix(IN_X, "OMP2", "X");
+  checkmatrix(IN_DtX, "OMP2", "DtX");
+  checkmatrix(IN_XtX, "OMP2", "XtX");
+  checkmatrix(IN_G, "OMP2", "G");
+  
+  checkscalar(IN_EPS, "OMP2", "EPSILON");
+  checkscalar(IN_SPARSE_G, "OMP2", "sparse_g");
+  checkscalar(IN_MSGDELTA, "OMP2", "msgdelta");
+  checkscalar(IN_MAXATOMS, "OMP2", "maxatoms");
+  checkscalar(IN_PROFILE, "OMP2", "profile");
+  
+  
+  /* get parameters */
+  
+  x = D = DtX = XtX = G = 0;
+  
+  if (!mxIsEmpty(IN_D))
+    D = mxGetPr(IN_D);
+  
+  if (!mxIsEmpty(IN_X))
+    x = mxGetPr(IN_X);
+  
+  if (!mxIsEmpty(IN_DtX))
+    DtX = mxGetPr(IN_DtX);
+  
+  if (!mxIsEmpty(IN_XtX))
+    XtX = mxGetPr(IN_XtX);
+  
+  if (!mxIsEmpty(IN_G))
+    G = mxGetPr(IN_G);
+  
+  eps = mxGetScalar(IN_EPS);
+  if ((int)(mxGetScalar(IN_SPARSE_G)+1e-2)) {
+    gmode = SPARSE_GAMMA;
+  }
+  else {
+    gmode = FULL_GAMMA;
+  }
+  msgdelta = mxGetScalar(IN_MSGDELTA);
+  if (mxGetScalar(IN_MAXATOMS) < -1e-5) {
+    maxatoms = -1;
+  }
+  else {
+    maxatoms = (int)(mxGetScalar(IN_MAXATOMS)+1e-2);
+  }
+  profile = (int)(mxGetScalar(IN_PROFILE)+1e-2);
+  
+  
+  /* check sizes */
+  
+  if (D && x) {
+    n = mxGetM(IN_D);
+    m = mxGetN(IN_D);
+    L = mxGetN(IN_X);
+    
+    if (mxGetM(IN_X) != n) {
+      mexErrMsgTxt("D and X have incompatible sizes.");
+    }
+    
+    if (G) {
+      if (mxGetN(IN_G)!=mxGetM(IN_G)) {
+        mexErrMsgTxt("G must be a square matrix.");
+      }
+      if (mxGetN(IN_G) != m) {
+        mexErrMsgTxt("D and G have incompatible sizes.");
+      }
+    }
+  }
+  
+  else if (DtX && XtX) {
+    m = mxGetM(IN_DtX);
+    L = mxGetN(IN_DtX);
+    
+    /* set n to an arbitrary value that is at least the max possible number of selected atoms */
+    
+    if (maxatoms>0) {
+      n = maxatoms;
+    }
+    else {
+      n = m;
+    }
+    
+    if ( !(mxGetM(IN_XtX)==L && mxGetN(IN_XtX)==1) && !(mxGetM(IN_XtX)==1 && mxGetN(IN_XtX)==L) ) {
+      mexErrMsgTxt("DtX and XtX have incompatible sizes.");
+    }
+    
+    if (mxGetN(IN_G)!=mxGetM(IN_G)) {
+      mexErrMsgTxt("G must be a square matrix.");
+    }
+    if (mxGetN(IN_G) != m) {
+      mexErrMsgTxt("DtX and G have incompatible sizes.");
+    }
+  }
+  
+  else {
+    mexErrMsgTxt("Either D and X, or DtX and XtX, must be specified.");
+  }
+  
+  
+  /* Do OMP! */
+  
+  GAMMA_OUT = ompcoreGabor(D, x, DtX, XtX, G, n, m, L, maxatoms, eps, gmode, profile, msgdelta, 1);
+  
+  return;
+}
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/ksvd utils/ompbox utils/omp2mexGabor.m	Thu Jul 21 14:07:41 2011 +0100
@@ -0,0 +1,23 @@
+%This is the Matlab interface to the OMP2 MEX implementation.
+%The function is not for independent use, only through omp2.m.
+
+
+%OMP2MEX Matlab interface to the OMP2 MEX implementation.
+%  GAMMA = OMP2MEXGABOR(D,X,DtX,XtX,G,EPSILON,SPARSE_G,MSGDELTA,MAXATOMS,PROFILE)
+%  invokes the OMP2 MEX function according to the specified parameters. Not
+%  all the parameters are required. Those among D, X, DtX, XtX and G which
+%  are not specified should be passed as [].
+%
+%  EPSILON - the target error.
+%  SPARSE_G - returns a sparse GAMMA when nonzero, full GAMMA when zero.
+%  MSGDELTA - the delay in secs between messages. Zero means no messages.
+%  MAXATOMS - the max number of atoms per signal, negative for no max.
+%  PROFILE - nonzero means that profiling information should be printed.
+
+
+%  Ron Rubinstein
+%  Computer Science Department
+%  Technion, Haifa 32000 Israel
+%  ronrubin@cs
+%
+%  April 2009
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/ksvd utils/ompbox utils/ompGabor.m	Thu Jul 21 14:07:41 2011 +0100
@@ -0,0 +1,180 @@
+function gamma = omp(varargin)
+%OMP Sparsity-constrained Orthogonal Matching Pursuit.
+%  GAMMA = OMP(D,X,G,T) solves the optimization problem
+%
+%       min  |X - D*GAMMA|_2     s.t.  |GAMMA|_0 <= T
+%      gamma
+%
+%  for each of the signals in X, using Batch Orthogonal Matching Pursuit.
+%  Here, D is a dictionary with normalized columns, X is a matrix
+%  containing column signals, T is the # of non-zeros in each signal
+%  representation, and G is the Gramm matrix D'*D. The output GAMMA is a
+%  matrix containing the sparse representations as its columns. 
+%
+%  GAMMA = OMP(D,X,[],T) performs the same operation, but without the
+%  matrix G, using OMP-Cholesky. This call produces the same output as
+%  Batch-OMP, but is significantly slower. Using this syntax is only
+%  recommended when available memory is too small to store G.
+%
+%  GAMMA = OMP(DtX,G,T) is the fastest implementation of OMP, but also
+%  requires the most memory. Here, DtX stores the projections D'*X. In this
+%  case Batch-OMP is used, but without having to compute D'*X in advance,
+%  which slightly improves runtime. Note that in general, the call
+%
+%    GAMMA = OMP(D'*X,G,T);
+%
+%  will be faster than the call
+%
+%    GAMMA = OMP(D,X,G,T);
+%
+%  due to optimized matrix multiplications in Matlab. However, when the
+%  entire matrix D'*X cannot be stored in memory, one of the other two
+%  versions can be used. Both compute D'*X for just one signal at a time,
+%  and thus require much less memory.
+%
+%  GAMMA = OMP(...,PARAM1,VAL1,PARAM2,VAL2,...) specifies additional
+%  parameters for OMP. Available parameters are:
+%
+%    'gammamode' - Specifies the representation mode for GAMMA. Can be
+%                  either 'full' or 'sparse', corresponding to a full or
+%                  sparse matrix, respectively. By default, GAMMA is
+%                  returned as a sparse matrix.
+%    'messages'  - Specifies whether progress messages should be displayed.
+%                  When positive, this is the number of seconds between
+%                  status prints. When negative, indicates that no messages
+%                  should be displayed (this is the default).
+%    'checkdict' - Specifies whether dictionary normalization should be
+%                  verified. When set to 'on' (default) the dictionary
+%                  atoms are verified to be of unit L2-norm. Setting this
+%                  parameter to 'off' disables verification and accelerates
+%                  function performance. Note that an unnormalized
+%                  dictionary will produce invalid results.
+%    'profile'   - Can be either 'on' or 'off'. When 'on', profiling
+%                  information is displayed at the end of the funciton
+%                  execution.
+%
+%
+%  Summary of OMP versions:
+%
+%    version      |   speed     |   memory
+%  --------------------------------------------------
+%   OMP(DtX,G,T)  |  very fast  |  very large
+%   OMP(D,X,G,T)  |  fast       |  moderate
+%   OMP(D,X,[],T) |  very slow  |  small
+%  --------------------------------------------------
+%
+%
+%  References:
+%  [1] M. Elad, R. Rubinstein, and M. Zibulevsky, "Efficient Implementation
+%      of the K-SVD Algorithm using Batch Orthogonal Matching Pursuit",
+%      Technical Report - CS, Technion, April 2008.
+%
+%  See also OMP2.
+
+
+%  Ron Rubinstein
+%  Computer Science Department
+%  Technion, Haifa 32000 Israel
+%  ronrubin@cs
+%
+%  April 2009
+
+
+% default options
+
+sparse_gamma = 1;
+msgdelta = -1;
+checkdict = 1;
+profile = 0;
+
+
+% determine number of parameters
+
+paramnum = 1;
+while (paramnum<=nargin && ~ischar(varargin{paramnum}))
+  paramnum = paramnum+1;
+end
+paramnum = paramnum-1;
+
+
+% parse options
+
+for i = paramnum+1:2:length(varargin)
+  paramname = varargin{i};
+  paramval = varargin{i+1};
+
+  switch lower(paramname)
+
+    case 'gammamode'
+      if (strcmpi(paramval,'sparse'))
+        sparse_gamma = 1;
+      elseif (strcmpi(paramval,'full'))
+        sparse_gamma = 0;
+      else
+        error('Invalid GAMMA mode');
+      end
+      
+    case 'messages'
+      msgdelta = paramval;
+
+    case 'checkdict'
+      if (strcmpi(paramval,'on'))
+        checkdict = 1;
+      elseif (strcmpi(paramval,'off'))
+        checkdict = 0;
+      else
+        error('Invalid checkdict option');
+      end
+
+    case 'profile'
+      if (strcmpi(paramval,'on'))
+        profile = 1;
+      elseif (strcmpi(paramval,'off'))
+        profile = 0;
+      else
+        error('Invalid profile mode');
+      end
+
+    otherwise
+      error(['Unknown option: ' paramname]);
+  end
+  
+end
+
+
+% determine call type
+
+if (paramnum==3)
+  DtX = varargin{1};
+  G = varargin{2};
+  T = varargin{3};
+  D = [];
+  X = [];
+elseif (paramnum==4)
+  D = varargin{1};
+  X = varargin{2};
+  G = varargin{3};
+  T = varargin{4};
+  DtX = [];
+else
+  error('Invalid number of parameters');
+end
+
+
+% verify dictionary normalization
+
+if (checkdict)
+  if (isempty(G))
+    atomnorms = sum(D.*D);
+  else
+    atomnorms = diag(G);
+  end
+  if (any(abs(atomnorms-1) > 1e-2))
+    error('Dictionary columns must be normalized to unit length');
+  end
+end
+
+
+% omp
+
+gamma = ompmexGabor(D,X,DtX,G,T,sparse_gamma,msgdelta,profile);
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/ksvd utils/ompbox utils/ompcore.c	Thu Jul 21 14:07:41 2011 +0100
@@ -0,0 +1,409 @@
+/**************************************************************************
+ *
+ * File name: ompcore.c
+ *
+ * Ron Rubinstein
+ * Computer Science Department
+ * Technion, Haifa 32000 Israel
+ * ronrubin@cs
+ *
+ * Last Updated: 25.8.2009
+ *
+ *************************************************************************/
+
+
+#include "ompcore.h"
+#include "omputils.h"
+#include "ompprof.h"
+#include "myblas.h"
+#include <math.h>
+#include <string.h>
+
+
+
+/******************************************************************************
+ *                                                                            *
+ *                           Batch-OMP Implementation                         *
+ *                                                                            *
+ ******************************************************************************/  
+
+mxArray* ompcore(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, signum, pos, *ind, *gammaIr, *gammaJc, gamma_count;
+  mwSize allocated_coefs, allocated_cols;
+  int DtX_specified, XtX_specified, batchomp, standardomp, *selected_atoms;
+  double *alpha, *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 ***/
+  
+  alpha = (double*)mxMalloc(m*sizeof(double));        /* contains D'*residual */
+  ind = (mwIndex*)mxMalloc(n*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);
+      }
+      
+      
+      /* initialize alpha := DtX */
+      
+      memcpy(alpha, 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 */
+      
+      pos = maxabs(alpha, m);
+      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;
+      }
+      
+      
+      /* mark selected atom */
+      
+      ind[i] = pos;
+      selected_atoms[pos] = 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*m, m*sizeof(double));
+      }
+      else {
+        memcpy(Dsub+i*n, D+pos*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, 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, alpha, n, m);            /* compute alpha := 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(alpha, DtX + m*signum*DtX_specified, m*sizeof(double));    /* set alpha = D'*x */
+        vec_sum(-1, tempvec1, alpha, m);                                  /* compute alpha := alpha - 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(alpha);
+  
+  return Gamma;
+}
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/ksvd utils/ompbox utils/ompcore.h	Thu Jul 21 14:07:41 2011 +0100
@@ -0,0 +1,80 @@
+/**************************************************************************
+ *
+ * File name: ompcore.h
+ *
+ * Ron Rubinstein
+ * Computer Science Department
+ * Technion, Haifa 32000 Israel
+ * ronrubin@cs
+ *
+ * Last Updated: 18.8.2009
+ *
+ * Contains the core implementation of Batch-OMP / OMP-Cholesky.
+ *
+ *************************************************************************/
+
+
+#ifndef __OMP_CORE_H__
+#define __OMP_CORE_H__
+
+
+#include "mex.h"
+
+
+
+/**************************************************************************
+ * Perform Batch-OMP or OMP-Cholesky on a specified set of signals, using
+ * either a fixed number of atoms or an error bound.
+ *
+ * Parameters (not all required):
+ *
+ *   D - the dictionary, of size n X m
+ *   x - the signals, of size n X L
+ *   DtX - D'*x, of size m X L
+ *   XtX - squared norms of the signals in x, sum(x.*x), of length L
+ *   G - D'*D, of size m X m
+ *   T - target sparsity, or maximal number of atoms for error-based OMP
+ *   eps - target residual norm for error-based OMP
+ *   gamma_mode - one of the constants FULL_GAMMA or SPARSE_GAMMA
+ *   profile - if non-zero, profiling info is printed
+ *   msg_delta - positive: the # of seconds between status prints, otherwise: nothing is printed
+ *   erroromp - if nonzero indicates error-based OMP, otherwise fixed sparsity OMP
+ *
+ * Usage:
+ *
+ *   The function can be called using different parameters, and will have
+ *   different complexity depending on the parameters specified. Arrays which
+ *   are not specified should be passed as null (0). When G is specified, 
+ *   Batch-OMP is performed. Otherwise, OMP-Cholesky is performed.
+ *
+ *   Fixed-sparsity usage:
+ *   ---------------------
+ *   Either DtX, or D and x, must be specified. Specifying DtX is more efficient.
+ *   XtX does not need to be specified.
+ *   When D and x are specified, G is not required. However, not providing G
+ *   will significantly degrade efficiency.
+ *   The number of atoms must be specified in T. The value of eps is ignored.
+ *   Finally, set erroromp to 0.
+ *
+ *   Error-OMP usage:
+ *   ----------------
+ *   Either DtX and Xtx, or D and x, must be specified. Specifying DtX and XtX
+ *   is more efficient.
+ *   When D and x are specified, G is not required. However, not providing G
+ *   will significantly degrade efficiency.
+ *   The target error must be specified in eps. A hard limit on the number
+ *   of atoms can also be specified via the parameter T. Otherwise, T should 
+ *   be negative. Finally, set erroromp to nonzero.
+ *
+ *
+ * Returns: 
+ *   An mxArray containing the sparse representations of the signals in x
+ *   (allocated using the appropriate mxCreateXXX() function).
+ *   The array is either full or sparse, depending on gamma_mode.
+ *
+ **************************************************************************/
+mxArray* ompcore(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);
+
+
+#endif
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/ksvd utils/ompbox utils/ompcoreGabor.c	Thu Jul 21 14:07:41 2011 +0100
@@ -0,0 +1,462 @@
+/**************************************************************************
+ * 
+ * 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++) {
+    for (j=0; j<n; j++) {
+          D1D2[i] += D1[i*n+j];
+    }
+    n12[i]=1/(1-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++){
+          for (j=0; j<n; j++){
+                  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]);
+          }
+      }
+      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;
+}
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/ksvd utils/ompbox utils/ompcoreGabor.h	Thu Jul 21 14:07:41 2011 +0100
@@ -0,0 +1,80 @@
+/**************************************************************************
+ *
+ * File name: ompcore.h
+ *
+ * Ron Rubinstein
+ * Computer Science Department
+ * Technion, Haifa 32000 Israel
+ * ronrubin@cs
+ *
+ * Last Updated: 18.8.2009
+ *
+ * Contains the core implementation of Batch-OMP / OMP-Cholesky.
+ *
+ *************************************************************************/
+
+
+#ifndef __OMP_CORE_H__
+#define __OMP_CORE_H__
+
+
+#include "mex.h"
+
+
+
+/**************************************************************************
+ * Perform Batch-OMP or OMP-Cholesky on a specified set of signals, using
+ * either a fixed number of atoms or an error bound.
+ *
+ * Parameters (not all required):
+ *
+ *   D - the dictionary, of size n X m
+ *   x - the signals, of size n X L
+ *   DtX - D'*x, of size m X L
+ *   XtX - squared norms of the signals in x, sum(x.*x), of length L
+ *   G - D'*D, of size m X m
+ *   T - target sparsity, or maximal number of atoms for error-based OMP
+ *   eps - target residual norm for error-based OMP
+ *   gamma_mode - one of the constants FULL_GAMMA or SPARSE_GAMMA
+ *   profile - if non-zero, profiling info is printed
+ *   msg_delta - positive: the # of seconds between status prints, otherwise: nothing is printed
+ *   erroromp - if nonzero indicates error-based OMP, otherwise fixed sparsity OMP
+ *
+ * Usage:
+ *
+ *   The function can be called using different parameters, and will have
+ *   different complexity depending on the parameters specified. Arrays which
+ *   are not specified should be passed as null (0). When G is specified, 
+ *   Batch-OMP is performed. Otherwise, OMP-Cholesky is performed.
+ *
+ *   Fixed-sparsity usage:
+ *   ---------------------
+ *   Either DtX, or D and x, must be specified. Specifying DtX is more efficient.
+ *   XtX does not need to be specified.
+ *   When D and x are specified, G is not required. However, not providing G
+ *   will significantly degrade efficiency.
+ *   The number of atoms must be specified in T. The value of eps is ignored.
+ *   Finally, set erroromp to 0.
+ *
+ *   Error-OMP usage:
+ *   ----------------
+ *   Either DtX and Xtx, or D and x, must be specified. Specifying DtX and XtX
+ *   is more efficient.
+ *   When D and x are specified, G is not required. However, not providing G
+ *   will significantly degrade efficiency.
+ *   The target error must be specified in eps. A hard limit on the number
+ *   of atoms can also be specified via the parameter T. Otherwise, T should 
+ *   be negative. Finally, set erroromp to nonzero.
+ *
+ *
+ * Returns: 
+ *   An mxArray containing the sparse representations of the signals in x
+ *   (allocated using the appropriate mxCreateXXX() function).
+ *   The array is either full or sparse, depending on gamma_mode.
+ *
+ **************************************************************************/
+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);
+
+
+#endif
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/ksvd utils/ompbox utils/ompmex.c	Thu Jul 21 14:07:41 2011 +0100
@@ -0,0 +1,133 @@
+/**************************************************************************
+ *
+ * File name: ompmex.c
+ *
+ * Ron Rubinstein
+ * Computer Science Department
+ * Technion, Haifa 32000 Israel
+ * ronrubin@cs
+ *
+ * Last Updated: 18.8.2009
+ *
+ *************************************************************************/
+
+#include "ompcore.h"
+#include "omputils.h"
+#include "mexutils.h"
+
+
+/* Input Arguments */
+
+#define IN_D          prhs[0]
+#define IN_X          prhs[1]
+#define IN_DtX        prhs[2]
+#define IN_G          prhs[3]
+#define IN_T          prhs[4]
+#define IN_SPARSE_G   prhs[5]
+#define IN_MSGDELTA   prhs[6]
+#define IN_PROFILE    prhs[7]
+
+
+/* Output Arguments */
+
+#define	GAMMA_OUT     plhs[0]
+
+
+/***************************************************************************************/
+
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray*prhs[])
+
+{
+  double *D, *x, *DtX, *G, msgdelta;
+  int gmode, profile, T;
+  mwSize m, n, L;   /* D is n x m , X is n x L, DtX is m x L */
+  
+  
+  /* check parameters */
+  
+  checkmatrix(IN_D, "OMP", "D");
+  checkmatrix(IN_X, "OMP", "X");
+  checkmatrix(IN_DtX, "OMP", "DtX");
+  checkmatrix(IN_G, "OMP", "G");
+  
+  checkscalar(IN_T, "OMP", "T");
+  checkscalar(IN_SPARSE_G, "OMP", "sparse_g");
+  checkscalar(IN_MSGDELTA, "OMP", "msgdelta");
+  checkscalar(IN_PROFILE, "OMP", "profile");
+
+  
+  /* get parameters */
+  
+  x = D = DtX = G = 0;
+  
+  if (!mxIsEmpty(IN_D))
+    D = mxGetPr(IN_D);
+  
+  if (!mxIsEmpty(IN_X))
+    x = mxGetPr(IN_X);
+  
+  if (!mxIsEmpty(IN_DtX))
+    DtX = mxGetPr(IN_DtX);
+  
+  if (!mxIsEmpty(IN_G))
+    G = mxGetPr(IN_G);
+  
+  T = (int)(mxGetScalar(IN_T)+1e-2);
+  if ((int)(mxGetScalar(IN_SPARSE_G)+1e-2)) {
+    gmode = SPARSE_GAMMA;
+  }
+  else {
+    gmode = FULL_GAMMA;
+  }
+  msgdelta = mxGetScalar(IN_MSGDELTA);
+  profile = (int)(mxGetScalar(IN_PROFILE)+1e-2);
+  
+  
+  /* check sizes */
+  
+  if (D && x) {
+    n = mxGetM(IN_D);
+    m = mxGetN(IN_D);
+    L = mxGetN(IN_X);
+    
+    if (mxGetM(IN_X) != n) {
+      mexErrMsgTxt("D and X have incompatible sizes.");
+    }
+    
+    if (G) {
+      if (mxGetN(IN_G)!=mxGetM(IN_G)) {
+        mexErrMsgTxt("G must be a square matrix.");
+      }
+      if (mxGetN(IN_G) != m) {
+        mexErrMsgTxt("D and G have incompatible sizes.");
+      }
+    }
+  }
+  
+  else if (DtX) {
+    m = mxGetM(IN_DtX);
+    L = mxGetN(IN_DtX);
+    
+    n = T;  /* arbitrary - it is enough to assume signal length is T */
+    
+    if (mxGetN(IN_G)!=mxGetM(IN_G)) {
+      mexErrMsgTxt("G must be a square matrix.");
+    }
+    if (mxGetN(IN_G) != m) {
+      mexErrMsgTxt("DtX and G have incompatible sizes.");
+    }
+  }
+  
+  else {
+    mexErrMsgTxt("Either D and X, or DtX, must be specified.");
+  }
+  
+  
+  /* Do OMP! */
+  
+  GAMMA_OUT = ompcore(D, x, DtX, 0, G, n, m, L, T, 0, gmode, profile, msgdelta, 0);
+  
+  return;
+}
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/ksvd utils/ompbox utils/ompmex.m	Thu Jul 21 14:07:41 2011 +0100
@@ -0,0 +1,22 @@
+%This is the Matlab interface to the OMP MEX implementation.
+%The function is not for independent use, only through omp.m.
+
+
+%OMPMEX Matlab interface to the OMP MEX implementation.
+%  GAMMA = OMPMEX(D,X,DtX,G,L,SPARSE_G,MSGDELTA,PROFILE) invokes the OMP
+%  MEX function according to the specified parameters. Not all the
+%  parameters are required. Those among D, X, DtX and G which are not
+%  specified should be passed as [].
+%
+%  L - the target sparsity.
+%  SPARSE_G - returns a sparse GAMMA when nonzero, full GAMMA when zero.
+%  MSGDELTA - the delay in secs between messages. Zero means no messages.
+%  PROFILE - nonzero means that profiling information should be printed.
+
+
+%  Ron Rubinstein
+%  Computer Science Department
+%  Technion, Haifa 32000 Israel
+%  ronrubin@cs
+%
+%  April 2009
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/ksvd utils/ompbox utils/ompmexGabor.c	Thu Jul 21 14:07:41 2011 +0100
@@ -0,0 +1,133 @@
+/**************************************************************************
+ *
+ * File name: ompmex.c
+ *
+ * Ron Rubinstein
+ * Computer Science Department
+ * Technion, Haifa 32000 Israel
+ * ronrubin@cs
+ *
+ * Last Updated: 18.8.2009
+ *
+ *************************************************************************/
+
+#include "ompcoreGabor.h"
+#include "omputils.h"
+#include "mexutils.h"
+
+
+/* Input Arguments */
+
+#define IN_D          prhs[0]
+#define IN_X          prhs[1]
+#define IN_DtX        prhs[2]
+#define IN_G          prhs[3]
+#define IN_T          prhs[4]
+#define IN_SPARSE_G   prhs[5]
+#define IN_MSGDELTA   prhs[6]
+#define IN_PROFILE    prhs[7]
+
+
+/* Output Arguments */
+
+#define	GAMMA_OUT     plhs[0]
+
+
+/***************************************************************************************/
+
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray*prhs[])
+
+{
+  double *D, *x, *DtX, *G, msgdelta;
+  int gmode, profile, T;
+  mwSize m, n, L;   /* D is n x m , X is n x L, DtX is m x L */
+  
+  
+  /* check parameters */
+  
+  checkmatrix(IN_D, "OMP", "D");
+  checkmatrix(IN_X, "OMP", "X");
+  checkmatrix(IN_DtX, "OMP", "DtX");
+  checkmatrix(IN_G, "OMP", "G");
+  
+  checkscalar(IN_T, "OMP", "T");
+  checkscalar(IN_SPARSE_G, "OMP", "sparse_g");
+  checkscalar(IN_MSGDELTA, "OMP", "msgdelta");
+  checkscalar(IN_PROFILE, "OMP", "profile");
+
+  
+  /* get parameters */
+  
+  x = D = DtX = G = 0;
+  
+  if (!mxIsEmpty(IN_D))
+    D = mxGetPr(IN_D);
+  
+  if (!mxIsEmpty(IN_X))
+    x = mxGetPr(IN_X);
+  
+  if (!mxIsEmpty(IN_DtX))
+    DtX = mxGetPr(IN_DtX);
+  
+  if (!mxIsEmpty(IN_G))
+    G = mxGetPr(IN_G);
+  
+  T = (int)(mxGetScalar(IN_T)+1e-2);
+  if ((int)(mxGetScalar(IN_SPARSE_G)+1e-2)) {
+    gmode = SPARSE_GAMMA;
+  }
+  else {
+    gmode = FULL_GAMMA;
+  }
+  msgdelta = mxGetScalar(IN_MSGDELTA);
+  profile = (int)(mxGetScalar(IN_PROFILE)+1e-2);
+  
+  
+  /* check sizes */
+  
+  if (D && x) {
+    n = mxGetM(IN_D);
+    m = mxGetN(IN_D);
+    L = mxGetN(IN_X);
+    
+    if (mxGetM(IN_X) != n) {
+      mexErrMsgTxt("D and X have incompatible sizes.");
+    }
+    
+    if (G) {
+      if (mxGetN(IN_G)!=mxGetM(IN_G)) {
+        mexErrMsgTxt("G must be a square matrix.");
+      }
+      if (mxGetN(IN_G) != m) {
+        mexErrMsgTxt("D and G have incompatible sizes.");
+      }
+    }
+  }
+  
+  else if (DtX) {
+    m = mxGetM(IN_DtX);
+    L = mxGetN(IN_DtX);
+    
+    n = T;  /* arbitrary - it is enough to assume signal length is T */
+    
+    if (mxGetN(IN_G)!=mxGetM(IN_G)) {
+      mexErrMsgTxt("G must be a square matrix.");
+    }
+    if (mxGetN(IN_G) != m) {
+      mexErrMsgTxt("DtX and G have incompatible sizes.");
+    }
+  }
+  
+  else {
+    mexErrMsgTxt("Either D and X, or DtX, must be specified.");
+  }
+  
+  
+  /* Do OMP! */
+  
+  GAMMA_OUT = ompcoreGabor(D, x, DtX, 0, G, n, m, L, T, 0, gmode, profile, msgdelta, 0);
+  
+  return;
+}
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/ksvd utils/ompbox utils/ompmexGabor.m	Thu Jul 21 14:07:41 2011 +0100
@@ -0,0 +1,22 @@
+%This is the Matlab interface to the OMP MEX implementation.
+%The function is not for independent use, only through omp.m.
+
+
+%OMPMEX Matlab interface to the OMP MEX implementation.
+%  GAMMA = OMPMEXGABOR(D,X,DtX,G,L,SPARSE_G,MSGDELTA,PROFILE) invokes the OMP
+%  MEX function according to the specified parameters. Not all the
+%  parameters are required. Those among D, X, DtX and G which are not
+%  specified should be passed as [].
+%
+%  L - the target sparsity.
+%  SPARSE_G - returns a sparse GAMMA when nonzero, full GAMMA when zero.
+%  MSGDELTA - the delay in secs between messages. Zero means no messages.
+%  PROFILE - nonzero means that profiling information should be printed.
+
+
+%  Ron Rubinstein
+%  Computer Science Department
+%  Technion, Haifa 32000 Israel
+%  ronrubin@cs
+%
+%  April 2009
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/ksvd utils/ompbox utils/ompprof.c	Thu Jul 21 14:07:41 2011 +0100
@@ -0,0 +1,113 @@
+/**************************************************************************
+ *
+ * File name: ompprof.c
+ *
+ * Ron Rubinstein
+ * Computer Science Department
+ * Technion, Haifa 32000 Israel
+ * ronrubin@cs
+ *
+ * Last Updated: 11.4.2009
+ *
+ *************************************************************************/
+
+
+#include "ompprof.h"
+
+
+/* initialize profiling information */
+
+void initprofdata(profdata *pd)
+{
+  pd->DtX_time = 0;
+  pd->XtX_time = 0;
+  pd->DtR_time = 0;
+  pd->maxabs_time = 0;
+  pd->DtD_time = 0;
+  pd->Lchol_time = 0;
+  pd->compcoef_time = 0;
+  pd->update_DtR_time = 0;
+  pd->update_resnorm_time = 0;
+  pd->compres_time = 0;
+  pd->indexsort_time = 0;
+  
+  pd->DtX_time_counted = 0;
+  pd->XtX_time_counted = 0;
+  pd->DtR_time_counted = 0;
+  pd->DtD_time_counted = 0;
+  pd->update_DtR_time_counted = 0;
+  pd->resnorm_time_counted = 0;
+  pd->compres_time_counted = 0;
+  pd->indexsort_time_counted = 0;
+  
+  pd->prevtime = clock();
+}
+
+
+/* add elapsed time to profiling data according to specified computation */
+
+void addproftime(profdata *pd, int comptype)
+{
+  switch(comptype) {
+    case DtX_TIME:            pd->DtX_time            += clock()-pd->prevtime; pd->DtX_time_counted = 1; break;
+    case XtX_TIME:            pd->XtX_time            += clock()-pd->prevtime; pd->XtX_time_counted = 1; break;
+    case DtR_TIME:            pd->DtR_time            += clock()-pd->prevtime; pd->DtR_time_counted = 1; break;
+    case DtD_TIME:            pd->DtD_time            += clock()-pd->prevtime; pd->DtD_time_counted = 1; break;
+    case COMPRES_TIME:        pd->compres_time        += clock()-pd->prevtime; pd->compres_time_counted = 1; break;
+    case UPDATE_DtR_TIME:     pd->update_DtR_time     += clock()-pd->prevtime; pd->update_DtR_time_counted = 1; break;
+    case UPDATE_RESNORM_TIME: pd->update_resnorm_time += clock()-pd->prevtime; pd->resnorm_time_counted = 1; break;
+    case INDEXSORT_TIME:      pd->indexsort_time      += clock()-pd->prevtime; pd->indexsort_time_counted = 1; break;
+    case MAXABS_TIME:         pd->maxabs_time         += clock()-pd->prevtime; break;
+    case LCHOL_TIME:          pd->Lchol_time          += clock()-pd->prevtime; break;
+    case COMPCOEF_TIME:       pd->compcoef_time       += clock()-pd->prevtime; break;
+  }
+  pd->prevtime = clock();
+}
+
+
+/* print profiling info */
+
+void printprofinfo(profdata *pd, int erroromp, int batchomp, int signum)
+{
+  clock_t tottime;
+  
+  tottime = pd->DtX_time + pd->XtX_time + pd->DtR_time + pd->DtD_time + pd->compres_time + pd->maxabs_time + 
+            pd->Lchol_time + pd->compcoef_time + pd->update_DtR_time + pd->update_resnorm_time + pd->indexsort_time;
+  
+  mexPrintf("\n\n*****  Profiling information for %s  *****\n\n", erroromp? "OMP2" : "OMP");
+  
+  mexPrintf("OMP mode: %s\n\n", batchomp? "Batch-OMP" : "OMP-Cholesky");
+  
+  mexPrintf("Total signals processed: %d\n\n", signum);
+  
+  if (pd->DtX_time_counted) {
+    mexPrintf("Compute DtX time:      %7.3lf seconds\n", pd->DtX_time/(double)CLOCKS_PER_SEC);
+  }
+  if (pd->XtX_time_counted) {
+    mexPrintf("Compute XtX time:      %7.3lf seconds\n", pd->XtX_time/(double)CLOCKS_PER_SEC);
+  }
+  mexPrintf("Max abs time:          %7.3lf seconds\n", pd->maxabs_time/(double)CLOCKS_PER_SEC);
+  if (pd->DtD_time_counted) {
+    mexPrintf("Compute DtD time:      %7.3lf seconds\n", pd->DtD_time/(double)CLOCKS_PER_SEC);
+  }
+  mexPrintf("Lchol update time:     %7.3lf seconds\n", pd->Lchol_time/(double)CLOCKS_PER_SEC);
+  mexPrintf("Compute coef time:     %7.3lf seconds\n", pd->compcoef_time/(double)CLOCKS_PER_SEC);
+  if (pd->compres_time_counted) {
+    mexPrintf("Compute R time:        %7.3lf seconds\n", pd->compres_time/(double)CLOCKS_PER_SEC);
+  }
+  if (pd->DtR_time_counted) {
+    mexPrintf("Compute DtR time:      %7.3lf seconds\n", pd->DtR_time/(double)CLOCKS_PER_SEC);
+  }
+  if (pd->update_DtR_time_counted) {
+    mexPrintf("Update DtR time:       %7.3lf seconds\n", pd->update_DtR_time/(double)CLOCKS_PER_SEC);
+  }
+  if (pd->resnorm_time_counted) {
+    mexPrintf("Update resnorm time:   %7.3lf seconds\n", pd->update_resnorm_time/(double)CLOCKS_PER_SEC);
+  }
+  if (pd->indexsort_time_counted) {
+    mexPrintf("Index sort time:       %7.3lf seconds\n", pd->indexsort_time/(double)CLOCKS_PER_SEC);
+  }
+  mexPrintf("---------------------------------------\n");
+  mexPrintf("Total time:            %7.3lf seconds\n\n", tottime/(double)CLOCKS_PER_SEC);
+}
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/ksvd utils/ompbox utils/ompprof.h	Thu Jul 21 14:07:41 2011 +0100
@@ -0,0 +1,106 @@
+/**************************************************************************
+ *
+ * File name: ompprof.h
+ *
+ * Ron Rubinstein
+ * Computer Science Department
+ * Technion, Haifa 32000 Israel
+ * ronrubin@cs
+ *
+ * Last Updated: 18.8.2009
+ *
+ * Collection of definitions and functions for profiling the OMP method.
+ *
+ *************************************************************************/
+
+
+#ifndef __OMP_PROF_H__
+#define __OMP_PROF_H__
+
+#include "mex.h"
+#include <time.h>
+
+
+
+/**************************************************************************
+ *
+ * Constants and data types.
+ *
+ **************************************************************************/
+
+
+/* constants denoting the various parts of the algorithm */
+
+enum { DtX_TIME, XtX_TIME, DtR_TIME, MAXABS_TIME, DtD_TIME, LCHOL_TIME, COMPCOEF_TIME, 
+       UPDATE_DtR_TIME, UPDATE_RESNORM_TIME, COMPRES_TIME, INDEXSORT_TIME };
+
+       
+       
+/* profiling data container with counters for each part of the algorithm */
+       
+typedef struct profdata 
+{
+  clock_t prevtime;  /* the time when last initialization/call to addproftime() was performed */
+  
+  clock_t DtX_time;
+  clock_t XtX_time;
+  clock_t DtR_time;
+  clock_t maxabs_time;
+  clock_t DtD_time;
+  clock_t Lchol_time;
+  clock_t compcoef_time;
+  clock_t update_DtR_time;
+  clock_t update_resnorm_time;
+  clock_t compres_time;
+  clock_t indexsort_time;
+  
+  /* flags indicating whether profiling data was gathered */
+  int DtX_time_counted;
+  int XtX_time_counted;
+  int DtR_time_counted;
+  int DtD_time_counted;
+  int update_DtR_time_counted;
+  int resnorm_time_counted;
+  int compres_time_counted;
+  int indexsort_time_counted;
+  
+} profdata;
+
+
+
+/**************************************************************************
+ *
+ * Initialize a profdata structure, zero all counters, and start its timer.
+ *
+ **************************************************************************/
+void initprofdata(profdata *pd);
+
+
+/**************************************************************************
+ *
+ * Add elapsed time from last call to addproftime(), or from initialization
+ * of profdata, to the counter specified by comptype. comptype must be one
+ * of the constants in the enumeration above.
+ *
+ **************************************************************************/
+void addproftime(profdata *pd, int comptype);
+
+
+/**************************************************************************
+ *
+ * Print the current contents of the counters in profdata.
+ *
+ * Parameters:
+ *   pd - the profdata to print
+ *   erroromp - indicates whether error-based (nonzero) or sparsity-based (zero)
+ *              omp was performed.
+ *   batchomp - indicates whether batch-omp (nonzero) or omp-cholesky (zero)
+ *              omp was performed.
+ *   signum   - number of signals processed by omp
+ *
+ **************************************************************************/
+void printprofinfo(profdata *pd, int erroromp, int batchomp, int signum);
+
+
+#endif
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/ksvd utils/ompbox utils/omputils.c	Thu Jul 21 14:07:41 2011 +0100
@@ -0,0 +1,89 @@
+/**************************************************************************
+ *
+ * File name: omputils.c
+ *
+ * Ron Rubinstein
+ * Computer Science Department
+ * Technion, Haifa 32000 Israel
+ * ronrubin@cs
+ *
+ * Last Updated: 18.8.2009
+ *
+ *************************************************************************/
+
+#include "omputils.h"
+#include <math.h>
+
+
+const char FULL_GAMMA_STR[] = "full";
+const char SPARSE_GAMMA_STR[] = "sparse";
+
+
+/* convert seconds to hours, minutes and seconds */
+
+void secs2hms(double sectot, int *hrs, int *mins, double *secs)
+{
+  *hrs = (int)(floor(sectot/3600)+1e-2);
+  sectot = sectot - 3600*(*hrs);
+  *mins = (int)(floor(sectot/60)+1e-2);
+  *secs = sectot - 60*(*mins);
+}
+
+
+/* quicksort, public-domain C implementation by Darel Rex Finley. */
+/* modification: sorts the array data[] as well, according to the values in the array vals[] */
+
+#define  MAX_LEVELS  300
+
+void quicksort(mwIndex vals[], double data[], mwIndex n) {
+  
+  long piv, beg[MAX_LEVELS], end[MAX_LEVELS], i=0, L, R, swap ;
+  double datapiv;
+  
+  beg[0]=0;
+  end[0]=n;
+  
+  while (i>=0) {
+    
+    L=beg[i]; 
+    R=end[i]-1;
+    
+    if (L<R) {
+      
+      piv=vals[L];
+      datapiv=data[L];
+      
+      while (L<R) 
+      {
+        while (vals[R]>=piv && L<R) 
+          R--;
+        if (L<R) {
+          vals[L]=vals[R];
+          data[L++]=data[R];
+        }
+        
+        while (vals[L]<=piv && L<R) 
+          L++;
+        if (L<R) { 
+          vals[R]=vals[L];
+          data[R--]=data[L];
+        }
+      }
+      
+      vals[L]=piv;
+      data[L]=datapiv;
+      
+      beg[i+1]=L+1;
+      end[i+1]=end[i];
+      end[i++]=L;
+      
+      if (end[i]-beg[i] > end[i-1]-beg[i-1]) {
+        swap=beg[i]; beg[i]=beg[i-1]; beg[i-1]=swap;
+        swap=end[i]; end[i]=end[i-1]; end[i-1]=swap;
+      }
+    }
+    else {
+      i--;
+    }
+  }
+}
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/ksvd utils/ompbox utils/omputils.h	Thu Jul 21 14:07:41 2011 +0100
@@ -0,0 +1,77 @@
+/**************************************************************************
+ *
+ * File name: omputils.h
+ *
+ * Ron Rubinstein
+ * Computer Science Department
+ * Technion, Haifa 32000 Israel
+ * ronrubin@cs
+ *
+ * Last Updated: 18.8.2009
+ *
+ * Utility definitions and functions for the OMP library.
+ *
+ *************************************************************************/
+
+
+#ifndef __OMP_UTILS_H__
+#define __OMP_UTILS_H__
+
+#include "mex.h"
+
+
+/* constants for the representation mode of gamma */
+
+extern const char FULL_GAMMA_STR[];      /* "full" */
+extern const char SPARSE_GAMMA_STR[];    /* "sparse" */
+
+
+#define FULL_GAMMA 1
+#define SPARSE_GAMMA 2
+#define INVALID_MODE 3
+
+
+
+/**************************************************************************
+ * Memory management for OMP2.
+ *
+ * GAMMA_INC_FACTOR:
+ * The matrix GAMMA is allocated with sqrt(n)/2 coefficients per signal,
+ * for a total of nzmax = L*sqrt(n)/2 nonzeros. Whenever GAMMA needs to be
+ * increased, it is increased by a factor of GAMMA_INC_FACTOR.
+ *
+ * MAT_INC_FACTOR:
+ * The matrices Lchol, Gsub and Dsub are allocated with sqrt(n)/2
+ * columns each. If additional columns are needed, this number is 
+ * increased by a factor of MAT_INC_FACTOR.
+ **************************************************************************/
+
+#define GAMMA_INC_FACTOR (1.4)
+#define MAT_INC_FACTOR (1.6)
+
+
+
+/**************************************************************************
+ * Convert number of seconds to hour, minute and second representation.
+ *
+ * Parameters:
+ *   sectot - total number of seconds
+ *   hrs, mins, secs - output hours (whole) and minutes (whole) and seconds
+ *
+ **************************************************************************/
+void secs2hms(double sectot, int *hrs, int *mins, double *secs);
+
+
+
+/**************************************************************************
+ * QuickSort - public-domain C implementation by Darel Rex Finley.
+ *
+ * Modified to sort both the array vals[] and the array data[] according 
+ * to the values in the array vals[].
+ *
+ **************************************************************************/
+void quicksort(mwIndex vals[], double data[], mwIndex n);
+
+
+#endif
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/ksvd utils/ompbox utils/printf.m	Thu Jul 21 14:07:41 2011 +0100
@@ -0,0 +1,26 @@
+function str = printf(varargin)
+%PRINTF Print formatted text to screen.
+%  PRINTF(FMT,VAL1,VAL2,...) formats the data in VAL1,VAL2,... according to
+%  the format string FMT, and prints the result to the screen.
+%
+%  The call to PRINTF(FMT,VAL1,VAL2,...) simply invokes the call
+%  DISP(SPRINTF(FMT,VAL1,VAL2,...)). For a complete description of the
+%  format string options see function SPRINTF.
+%
+%  STR = PRINTF(...) also returns the formatted string.
+
+
+%  Ron Rubinstein
+%  Computer Science Department
+%  Technion, Haifa 32000 Israel
+%  ronrubin@cs
+%
+%  April 2008
+
+
+if (nargout>0)
+  str = sprintf(varargin{:});
+  disp(str);
+else
+  disp(sprintf(varargin{:}));
+end