# HG changeset patch # User idamnjanovic # Date 1300121547 0 # Node ID 217a33ac374efc0f7b9f3fea114419d66b28d5ea # Parent d5155eaa3f6889aa2851190fcc8c58928fe0aaae diff -r d5155eaa3f68 -r 217a33ac374e Problems/private/genSampling.m --- /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 +% (mlustig@stanford.edu), 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)) + + +/* 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 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 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=1; --i) { + rhs = b[i-1]; + for (j=i; j=1; --i) { + rhs = b[i-1]; + for (j=i; j + + + +/************************************************************************** + * 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 + diff -r d5155eaa3f68 -r 217a33ac374e Problems/private/omp2mex.c --- /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; +} diff -r d5155eaa3f68 -r 217a33ac374e Problems/private/omp2mex.m --- /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 diff -r d5155eaa3f68 -r 217a33ac374e Problems/private/ompcore.h --- /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 diff -r d5155eaa3f68 -r 217a33ac374e Problems/private/ompmex.c --- /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; +} + diff -r d5155eaa3f68 -r 217a33ac374e Problems/private/ompmex.m --- /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 diff -r d5155eaa3f68 -r 217a33ac374e Problems/private/ompprof.c --- /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); +} + diff -r d5155eaa3f68 -r 217a33ac374e Problems/private/ompprof.h --- /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 + + + +/************************************************************************** + * + * 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 + diff -r d5155eaa3f68 -r 217a33ac374e Problems/private/omputils.c --- /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 + + +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=piv && L 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--; + } + } +} diff -r d5155eaa3f68 -r 217a33ac374e Problems/private/omputils.h --- /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 + diff -r d5155eaa3f68 -r 217a33ac374e Problems/private/pwsmoothfield.m --- /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 +% (swright@cs.wisc.edu), 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 + diff -r d5155eaa3f68 -r 217a33ac374e Problems/private/scalarToRGB.m --- /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 diff -r d5155eaa3f68 -r 217a33ac374e Problems/private/thumbPlot.m --- /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 diff -r d5155eaa3f68 -r 217a33ac374e Problems/private/thumbwrite.m --- /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); diff -r d5155eaa3f68 -r 217a33ac374e Problems/private/updateFigure.m --- /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