changeset 51:217a33ac374e

(none)

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Problems/private/genSampling.m	Mon Mar 14 16:52:27 2011 +0000
@@ -0,0 +1,53 @@
+function [minIntrVec,stat,actpctg] = genSampling(pdf,iter,tol)
+
+%[mask,stat,N] = genSampling(pdf,iter,tol)
+%
+% a monte-carlo algorithm to generate a sampling pattern with 
+% minimum peak interference. The number of samples will be
+% sum(pdf) +- tol
+%
+%	pdf - probability density function to choose samples from
+%	iter - number of tries
+%	tol  - the deviation from the desired number of samples in samples
+%
+% returns:
+%	mask - sampling pattern
+%	stat - vector of min interferences measured each try
+%	actpctg    - actual undersampling factor
+%
+%	(c) Michael Lustig 2007
+
+%   This file is used with the kind permission of Michael Lustig
+%   ([e-mail]), and originally appeared in the
+%   SparseMRI toolbox, http://www.stanford.edu/~mlustig/ .
+%
+% $Id: genSampling.m 1040 2008-06-26 20:29:02Z ewout78 $
+
+% h = waitbar(0);
+
+pdf(find(pdf>1)) = 1;
+K = sum(pdf(:));
+
+minIntr = 1e99;
+minIntrVec = zeros(size(pdf));
+
+for n=1:iter
+	tmp = zeros(size(pdf));
+	while abs(sum(tmp(:)) - K) > tol
+		tmp = rand(size(pdf))<pdf;
+	end
+	
+	TMP = ifft2(tmp./pdf);
+	if max(abs(TMP(2:end))) < minIntr
+		minIntr = max(abs(TMP(2:end)));
+		minIntrVec = tmp;
+	end
+	stat(n) = max(abs(TMP(2:end)));
+      % waitbar(n/iter,h);
+end
+
+actpctg = sum(minIntrVec(:))/prod(size(minIntrVec));
+
+% close(h);
+
+

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Problems/private/myblas.c	Mon Mar 14 16:52:27 2011 +0000
@@ -0,0 +1,663 @@
+/**************************************************************************
+ *
+ * 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*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/Problems/private/myblas.h	Mon Mar 14 16:52:27 2011 +0000
@@ -0,0 +1,492 @@
+/**************************************************************************
+ *
+ * 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);
+
+
+
+/**************************************************************************
+ * 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/Problems/private/omp2mex.c	Mon Mar 14 16:52:27 2011 +0000
@@ -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/Problems/private/omp2mex.m	Mon Mar 14 16:52:27 2011 +0000
@@ -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/Problems/private/ompcore.h	Mon Mar 14 16:52:27 2011 +0000
@@ -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/Problems/private/ompmex.c	Mon Mar 14 16:52:27 2011 +0000
@@ -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/Problems/private/ompmex.m	Mon Mar 14 16:52:27 2011 +0000
@@ -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/Problems/private/ompprof.c	Mon Mar 14 16:52:27 2011 +0000
@@ -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/Problems/private/ompprof.h	Mon Mar 14 16:52:27 2011 +0000
@@ -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/Problems/private/omputils.c	Mon Mar 14 16:52:27 2011 +0000
@@ -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/Problems/private/omputils.h	Mon Mar 14 16:52:27 2011 +0000
@@ -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/Problems/private/pwsmoothfield.m	Mon Mar 14 16:52:27 2011 +0000
@@ -0,0 +1,32 @@
+function f = pwsmoothfield(n,var,alpha)
+% PWSMOOTHFIELD(N,VAR,ALPHA)
+%   Generate an image of piecewise smooth filtered white noise
+%
+%   N = sidelength of the field in pixels
+%   VAR = variance of original white nose
+%   ALPHA = fraction of FFT coefficents to keep
+%
+%   Returns an N-by-N array
+%
+%   This file is used with the kind permission of Stephen J. Wright
+%   ([e-mail]), and was originally included in the GPSR
+%   v1.0 distribution: http://www.lx.it.pt/~mtf/GPSR .
+
+% $Id: pwsmoothfield.m 1040 2008-06-26 20:29:02Z ewout78 $
+
+f = sqrt(var)*randn(n);
+F = fft2(f);
+a = ceil(n*alpha/2);
+b = fix(n*(1-alpha));
+F(a+1:a+b,:) = 0;
+F(:,a+1:a+b) = 0;
+f = real(ifft2(F));
+
+for i = 1:n
+    for j = 1:n
+        if (j/n >= 15*(i/n - 0.5)^3 + 0.4)
+            f(i,j) = f(i,j) + sqrt(var);
+        end
+    end
+end
+

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Problems/private/scalarToRGB.m	Mon Mar 14 16:52:27 2011 +0000
@@ -0,0 +1,18 @@
+function s = scalarToRGB(x,colors)
+% input values are assumed to lie between 0 and 1
+
+%   Copyright 2008, Ewout van den Berg and Michael P. Friedlander
+%   http://www.cs.ubc.ca/labs/scl/sparco
+%   $Id: scalarToRGB.m 1040 2008-06-26 20:29:02Z ewout78 $
+
+l  = size(colors,1);
+m  = size(x,1);
+n  = size(x,2);
+s  = zeros(m,n,3);
+
+for i=1:m
+  for j=1:n
+     idx = max(1,min(l,1+floor((l-1) * x(i,j))));
+     s(i,j,:) = colors(idx,:);
+  end
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Problems/private/thumbPlot.m	Mon Mar 14 16:52:27 2011 +0000
@@ -0,0 +1,38 @@
+function P = thumbPlot(P,x,y,color)
+
+%   Copyright 2008, Ewout van den Berg and Michael P. Friedlander
+%   http://www.cs.ubc.ca/labs/scl/sparco
+%   $Id: thumbPlot.m 1040 2008-06-26 20:29:02Z ewout78 $
+
+m = size(P,1);
+n = size(P,2);
+if (size(P,3) == 0) & (length(color) == 3)
+  % Convert to gray-scale
+  color = 0.30*color(1) + 0.59*color(2) + 0.11*color(3);
+end
+
+mnx = min(x);  % Minimum x
+mxx = max(x);  % Maximum x
+mny = min(y);  % Minimum y
+mxy = max(y);  % Maximum y
+dy  = (mxy - mny) * 0.1;   % Offset on vertical axis
+sx  = (mxx - mnx) * 1.0;   % Scale of horizontal axis
+sy  = (mxy - mny) * 1.2;   % Scale of vertical axis
+
+if (sx < 1e-6), sx = 1; end
+if (sy < 1e-6), sy = 1; end
+
+for i=1:length(x)-1
+   x0 = floor(1 + (n-1) * (x(i  ) - mnx) / sx);
+   x1 = floor(1 + (n-1) * (x(i+1) - mnx) / sx);
+   y0 = floor(    (n-1) * (y(i  ) - mny + dy) / sy);
+   y1 = floor(    (n-1) * (y(i+1) - mny + dy) / sy);
+   
+   samples = 1+2*max(abs(x1-x0)+1,abs(y1-y0)+1);
+   c       = linspace(0,1,samples);
+   idx     = round((1-c)*x0 + c*x1);
+   idy     = n - round((1-c)*y0 + c*y1);
+   for j=1:samples
+      P(idy(j),idx(j),:) = color;
+   end
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Problems/private/thumbwrite.m	Mon Mar 14 16:52:27 2011 +0000
@@ -0,0 +1,18 @@
+function thumbwrite(data,name,opts)
+
+
+%   Copyright 2008, Ewout van den Berg and Michael P. Friedlander
+%   http://www.cs.ubc.ca/labs/scl/sparco
+%   $Id: thumbwrite.m 1040 2008-06-26 20:29:02Z ewout78 $
+
+%
+% data       Thumbnail data in range 0-1
+% name       Name of file (no extension)
+% opts
+%   .thumbtype  Type of image (png, eps, ps, ...)
+%   .thumbdir   Output directory
+%
+
+[type,ext] = getFigureExt(opts.thumbtype);
+data = round(data * 255) / 255;
+imwrite(data,[opts.thumbpath,name,'.',ext],type);

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Problems/private/updateFigure.m	Mon Mar 14 16:52:27 2011 +0000
@@ -0,0 +1,93 @@
+function updateFigure(opts, figTitle, figFilename)
+
+%   Copyright 2008, Ewout van den Berg and Michael P. Friedlander
+%   http://www.cs.ubc.ca/labs/scl/sparco
+%   $Id: updateFigure.m 1040 2008-06-26 20:29:02Z ewout78 $
+
+% Ensure default values are available
+opts.linewidth  = getOption(opts,'linewidth', []);
+opts.fontsize   = getOption(opts,'fontsize',  []);
+opts.markersize = getOption(opts,'markersize',[]);
+
+% Output the plots
+if opts.update
+  % Set the line width, font size and marker size
+  chld = [gca; get(gca,'Children')];
+  lnwd = ones(length(chld),1) * NaN;
+  fnts = ones(length(chld),1) * NaN;
+  mrks = ones(length(chld),1) * NaN;
+  for i=1:length(chld)
+    conf = get(chld(i));
+    if ~isempty(opts.linewidth) && isfield(conf,'LineWidth')
+      lnwd(i) = get(chld(i),'LineWidth');
+      if (lnwd(i) == 0.5) % Default
+        set(chld(i),'Linewidth',opts.linewidth);
+      end
+    end
+    if ~isempty(opts.fontsize) && isfield(conf,'FontSize')
+      fnts(i) = get(chld(i),'FontSize');
+      if (fnts(i) == 10) % Default
+        set(chld(i),'FontSize',opts.fontsize);
+      end
+    end
+    if ~isempty(opts.markersize) && isfield(conf,'MarkerSize')
+      mrks(i) = get(chld(i),'MarkerSize');
+      if (mrks(i) == 6) % Default
+        set(chld(i),'MarkerSize',opts.markersize);
+      end
+    end
+  end
+  
+  for i=1:length(opts.figtype)
+     updateFigureType(opts.update, 0, opts.figtype{i}, ...
+                      opts.figpath, figTitle, figFilename);
+  end
+
+  % Restore the line-widths, font size
+  for i=1:length(chld)
+    if ~isnan(lnwd(i))
+      set(chld(i),'LineWidth',lnwd(i));
+    end
+    if ~isnan(fnts(i))
+      set(chld(i),'FontSize',fnts(i));
+    end
+    if ~isnan(mrks(i))
+      set(chld(i),'MarkerSize',mrks(i));
+    end
+  end
+    
+end
+
+% Show the plot
+if opts.show
+  updateFigureType(0,opts.show,'','',figTitle,'');
+end
+
+
+
+function updateFigureType(update,show,figtype,figpath,figTitle,figFilename)
+filename = [figpath,figFilename];
+
+switch lower(figtype)
+ case {'pdf'}
+   cmdPostprocess = sprintf('!pdfcrop %s.pdf %s.pdf >& /dev/null', ...
+                            filename, filename);
+ otherwise
+   cmdPostprocess = [];
+end
+
+[figtype,figext] = getFigureExt(figtype);
+
+% Print the figure for output (without title)
+if update
+  evalc(sprintf('print -d%s %s.%s;', figtype, filename, figext));
+  
+  if ~isempty(cmdPostprocess)
+    eval(cmdPostprocess);
+  end
+end
+
+% Add title if needed
+if show
+  title(figTitle);
+end