Mercurial > hg > smallbox
changeset 59:23f9dd7b9d78
(none)
| author | idamnjanovic |
|---|---|
| date | Mon, 14 Mar 2011 17:25:38 +0000 |
| parents | 66fe897f1bf4 |
| children | ad36f80e2ccf |
| files | Problems/private/add_dc.m Problems/private/addtocols.c Problems/private/addtocols.m Problems/private/addtocols.mexa64 Problems/private/col2imstep.c Problems/private/col2imstep.m Problems/private/collincomb.c Problems/private/collincomb.m Problems/private/completeOps.m Problems/private/countcover.m Problems/private/diag_ids.m Problems/private/dictdist.m Problems/private/genPDF.m Problems/private/genSampling.m Problems/private/im2colstep.c Problems/private/im2colstep.m Problems/private/imnormalize.m Problems/private/iswhole.m Problems/private/make.m Problems/private/mexutils.c Problems/private/mexutils.h Problems/private/myblas.c Problems/private/myblas.h Problems/private/normcols.m Problems/private/omp2mex.c Problems/private/omp2mex.m Problems/private/ompcore.c Problems/private/ompcore.h Problems/private/ompmex.c Problems/private/ompmex.m Problems/private/ompprof.c Problems/private/ompprof.h Problems/private/omputils.c Problems/private/omputils.h Problems/private/printf.m Problems/private/private/make.m Problems/private/private/mexutils.c Problems/private/private/mexutils.h Problems/private/private/printf.m Problems/private/pwsmoothfield.m Problems/private/reggrid.m Problems/private/remove_dc.m Problems/private/rowlincomb.c Problems/private/rowlincomb.m Problems/private/sampgrid.m Problems/private/scalarToRGB.m Problems/private/secs2hms.m Problems/private/spdiag.m Problems/private/sprow.c Problems/private/sprow.m Problems/private/thumbFromOp.m Problems/private/thumbPlot.m Problems/private/thumbwrite.m Problems/private/timerclear.m Problems/private/timereta.m Problems/private/timerinit.m Problems/private/updateFigure.m |
| diffstat | 57 files changed, 0 insertions(+), 5015 deletions(-) [+] |
line wrap: on
line diff
--- a/Problems/private/add_dc.m Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,33 +0,0 @@ -function x = add_dc(y,dc,columns) -%ADD_DC Add DC channel to signals. -% X = ADD_DC(Y,DC) adds the specified DC value to the (possibly -% multi-dimensional) signal Y, returning the result as X. DC should be a -% scalar value. -% -% X = ADD_DC(Y,DC,'columns') where Y is a 2D matrix and DC is an array of -% length size(Y,2), treats the columns of Y as individual 1D signals, -% adding to each one the corresponding DC value from the DC array. X is -% the same size as Y and contains the resulting signals. -% -% See also REMOVE_DC. - -% Ron Rubinstein -% Computer Science Department -% Technion, Haifa 32000 Israel -% ronrubin@cs -% -% April 2009 - - -if (nargin==3 && strcmpi(columns,'columns')), columns = 1; -else columns = 0; -end - -if (columns) - x = addtocols(y,dc); -else - x = y + dc; -end - - -
--- a/Problems/private/addtocols.c Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,85 +0,0 @@ -/************************************************************************** - * - * File name: addtocols.c - * - * Ron Rubinstein - * Computer Science Department - * Technion, Haifa 32000 Israel - * ronrubin@cs - * - * Last Updated: 19.4.2009 - * - *************************************************************************/ - - -#include "mex.h" - - -/* Input Arguments */ - -#define X_IN prhs[0] -#define V_IN prhs[1] - - -/* Output Arguments */ - -#define Y_OUT plhs[0] - - -void mexFunction(int nlhs, mxArray *plhs[], - int nrhs, const mxArray*prhs[]) - -{ - double *x, *y, *v, *xend; - mwSize m,n,m1,n1; - mwIndex counter; - - - /* Check for proper number of arguments */ - - if (nrhs != 2) { - mexErrMsgTxt("Two input arguments required."); - } else if (nlhs > 1) { - mexErrMsgTxt("Too many output arguments."); - } - - - /* Check the the input dimensions */ - - m = mxGetM(X_IN); - n = mxGetN(X_IN); - if (!mxIsDouble(X_IN) || mxIsComplex(X_IN) || mxGetNumberOfDimensions(X_IN)>2) { - mexErrMsgTxt("ADDTOCOLS requires that X be a double matrix."); - } - m1 = mxGetM(V_IN); - n1 = mxGetN(V_IN); - if (!mxIsDouble(V_IN) || mxIsComplex(V_IN) || (m1!=1 && n1!=1)) { - mexErrMsgTxt("ADDTOCOLS requires that V be a double vector."); - } - if ((m1==1 && n1!=n) || (n1==1 && m1!=n)) { - mexErrMsgTxt("Error in ADDTOCOLS: dimensions of V and X must agree."); - } - - - /* Create a matrix for the return argument */ - Y_OUT = mxCreateDoubleMatrix(m, n, mxREAL); - - - /* Assign pointers to the various parameters */ - x = mxGetPr(X_IN); - v = mxGetPr(V_IN); - y = mxGetPr(Y_OUT); - - - /* Do the actual computation */ - - xend = x+(m*n); - counter = 0; - while (x<xend) { - (*y) = (*x) + (*v); - y++; x++; counter++; - if (counter==m) {v++; counter=0;} - } - - return; -}
--- a/Problems/private/addtocols.m Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,13 +0,0 @@ -%ADDTOCOLS Add values to the columns of a matrix. -% Y=ADDTOCOLS(X,V) adds to each column of the MxN matrix X the -% corresponding value from the N-element vector V. -% -% See also NORMCOLS. - - -% Ron Rubinstein -% Computer Science Department -% Technion, Haifa 32000 Israel -% ronrubin@cs -% -% June 2005
--- a/Problems/private/col2imstep.c Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,146 +0,0 @@ -/************************************************************************** - * - * File name: col2imstep.c - * - * Ron Rubinstein - * Computer Science Department - * Technion, Haifa 32000 Israel - * ronrubin@cs - * - * Last Updated: 31.8.2009 - * - *************************************************************************/ - - -#include "mex.h" - - -/* Input Arguments */ - -#define B_IN prhs[0] -#define N_IN prhs[1] -#define SZ_IN prhs[2] -#define S_IN prhs[3] - - -/* Output Arguments */ - -#define X_OUT plhs[0] - - -void mexFunction(int nlhs, mxArray *plhs[], - int nrhs, const mxArray*prhs[]) - -{ - double *x, *b, *s; - mwSize sz[3], stepsize[3], n[3], ndims; - mwIndex i, j, k, l, m, t, blocknum; - - - /* Check for proper number of arguments */ - - if (nrhs < 3 || nrhs > 4) { - mexErrMsgTxt("Invalid number of input arguments."); - } else if (nlhs > 1) { - mexErrMsgTxt("Too many output arguments."); - } - - - /* Check the the input dimensions */ - - if (!mxIsDouble(B_IN) || mxIsComplex(B_IN) || mxGetNumberOfDimensions(B_IN)>2) { - mexErrMsgTxt("B should be a double matrix."); - } - if (!mxIsDouble(N_IN) || mxIsComplex(N_IN) || mxGetNumberOfDimensions(N_IN)>2) { - mexErrMsgTxt("Invalid output matrix size."); - } - ndims = mxGetM(N_IN)*mxGetN(N_IN); - if (ndims<2 || ndims>3) { - mexErrMsgTxt("Output matrix can only be 2-D or 3-D."); - } - if (!mxIsDouble(SZ_IN) || mxIsComplex(SZ_IN) || mxGetNumberOfDimensions(SZ_IN)>2 || mxGetM(SZ_IN)*mxGetN(SZ_IN)!=ndims) { - mexErrMsgTxt("Invalid block size."); - } - if (nrhs == 4) { - if (!mxIsDouble(S_IN) || mxIsComplex(S_IN) || mxGetNumberOfDimensions(S_IN)>2 || mxGetM(S_IN)*mxGetN(S_IN)!=ndims) { - mexErrMsgTxt("Invalid step size."); - } - } - - - /* Get parameters */ - - s = mxGetPr(N_IN); - if (s[0]<1 || s[1]<1 || (ndims==3 && s[2]<1)) { - mexErrMsgTxt("Invalid output matrix size."); - } - n[0] = (mwSize)(s[0] + 0.01); - n[1] = (mwSize)(s[1] + 0.01); - n[2] = ndims==3 ? (mwSize)(s[2] + 0.01) : 1; - - s = mxGetPr(SZ_IN); - if (s[0]<1 || s[1]<1 || (ndims==3 && s[2]<1)) { - mexErrMsgTxt("Invalid block size."); - } - sz[0] = (mwSize)(s[0] + 0.01); - sz[1] = (mwSize)(s[1] + 0.01); - sz[2] = ndims==3 ? (mwSize)(s[2] + 0.01) : 1; - - if (nrhs == 4) { - s = mxGetPr(S_IN); - if (s[0]<1 || s[1]<1 || (ndims==3 && s[2]<1)) { - mexErrMsgTxt("Invalid step size."); - } - stepsize[0] = (mwSize)(s[0] + 0.01); - stepsize[1] = (mwSize)(s[1] + 0.01); - stepsize[2] = ndims==3 ? (mwSize)(s[2] + 0.01) : 1; - } - else { - stepsize[0] = stepsize[1] = stepsize[2] = 1; - } - - if (n[0]<sz[0] || n[1]<sz[1] || (ndims==3 && n[2]<sz[2])) { - mexErrMsgTxt("Block size too large."); - } - - if (mxGetN(B_IN) != ((n[0]-sz[0])/stepsize[0]+1)*((n[1]-sz[1])/stepsize[1]+1)*((n[2]-sz[2])/stepsize[2]+1)) { - mexErrMsgTxt("Invalid number of columns in B. Please use IM2COLSTEP to compute B."); - } - - - /* Create a matrix for the return argument */ - - X_OUT = mxCreateNumericArray(ndims, n, mxDOUBLE_CLASS, mxREAL); - - - /* Assign pointers */ - - b = mxGetPr(B_IN); - x = mxGetPr(X_OUT); - - - /* Do the actual computation */ - - blocknum = 0; - - /* iterate over all blocks */ - for (k=0; k<=n[2]-sz[2]; k+=stepsize[2]) { - for (j=0; j<=n[1]-sz[1]; j+=stepsize[1]) { - for (i=0; i<=n[0]-sz[0]; i+=stepsize[0]) { - - /* add single block */ - for (m=0; m<sz[2]; m++) { - for (l=0; l<sz[1]; l++) { - for (t=0; t<sz[0]; t++) { - (x+(k+m)*n[0]*n[1]+(j+l)*n[0]+i)[t] += (b + blocknum*sz[0]*sz[1]*sz[2] + m*sz[0]*sz[1] + l*sz[0])[t]; - } - } - } - blocknum++; - - } - } - } - - return; -}
--- a/Problems/private/col2imstep.m Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,25 +0,0 @@ -%COL2IMSTEP Rearrange matrix columns into blocks. -% A = COL2IMSTEP(B,[MM NN],[N1 N2]) rearranges the columns of B into -% sliding N1-by-N2 blocks producing the matrix A of size MM-by-NN. B is -% usually the result of calling IM2COLSTEP(...) with a stepsize of 1, or -% using Matlab's IM2COL(..,'sliding'). Overlapping blocks are summed in A. -% -% A = COL2IMSTEP(B,[MM NN],[N1 N2],[S1 S2]) arranges the blocks in A with -% a step size of (S1,S2) between them. The first block is at A(1:N1,1:N2), -% and the rest are at A((1:N1)+i*S1,(1:N2)+j*S2). Overlapping blocks are -% summed in A. Note that B is usually the result of calling -% IM2COLSTEP(...) with a stepsize of [S1 S2]. -% -% A = IM2COLSTEP(B,[MM NN KK],[N1 N2 N3],[S1 S2 S3]) generates a 3-D -% output matrix A. The step size [S1 S2 S3] may be ommitted, and defaults -% to [1 1 1]. -% -% See also IM2COLSTEP, IM2COL, COUNTCOVER. - - -% Ron Rubinstein -% Computer Science Department -% Technion, Haifa 32000 Israel -% ronrubin@cs -% -% August 2009
--- a/Problems/private/collincomb.c Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,165 +0,0 @@ -/************************************************************************** - * - * File name: collincomb.c - * - * Ron Rubinstein - * Computer Science Department - * Technion, Haifa 32000 Israel - * ronrubin@cs - * - * Last Updated: 21.5.2009 - * - *************************************************************************/ - - -#include "mex.h" - - -/* Input Arguments */ - -#define A_IN prhs[0] -#define ROWS_IN prhs[1] -#define COLS_IN1 prhs[1] -#define COLS_IN2 prhs[2] -#define X_IN1 prhs[2] -#define X_IN2 prhs[3] - - -/* Output Arguments */ - -#define Y_OUT plhs[0] - - -void mexFunction(int nlhs, mxArray *plhs[], -int nrhs, const mxArray*prhs[]) - -{ - double *A, *x, *y, *rows, *cols; - mwSize m,n,m1,n1,m2,n2,rownum,colnum; - mwIndex col_id,*row_ids,i,j; - int rownumspecified=0; - - - /* Check for proper number of arguments */ - - if (nrhs!=3 && nrhs!=4) { - mexErrMsgTxt("Invalid number of arguments."); - } else if (nlhs > 1) { - mexErrMsgTxt("Too many output arguments."); - } - - - /* Check the input dimensions */ - - m = mxGetM(A_IN); - n = mxGetN(A_IN); - if (!mxIsDouble(A_IN) || mxIsComplex(A_IN) || mxGetNumberOfDimensions(A_IN)>2) { - mexErrMsgTxt("COLLINCOMB requires that A be a double matrix."); - } - - if (nrhs==3) { - - m1 = mxGetM(COLS_IN1); - n1 = mxGetN(COLS_IN1); - if (!mxIsDouble(COLS_IN1) || mxIsComplex(COLS_IN1) || (m1!=1 && n1!=1)) { - mexErrMsgTxt("COLLINCOMB requires that COLS be an index vector of type double."); - } - colnum = (m1 > n1) ? m1 : n1; /* the number of columns in the linear combination */ - - m2 = mxGetM(X_IN1); - n2 = mxGetN(X_IN1); - if (!mxIsDouble(X_IN1) || mxIsComplex(X_IN1) || (m2!=1 && n2!=1)) { - mexErrMsgTxt("COLLINCOMB requires that X be a double vector."); - } - - if (m2!=colnum && n2!=colnum) { - mexErrMsgTxt("The length of X does not match the number of columns in COLS."); - } - - rows = 0; - Y_OUT = mxCreateDoubleMatrix(m, 1, mxREAL); - cols = mxGetPr(COLS_IN1); - x = mxGetPr(X_IN1); - } - else { - - m1 = mxGetM(ROWS_IN); - n1 = mxGetN(ROWS_IN); - if (!mxIsDouble(ROWS_IN) || mxIsComplex(ROWS_IN) || (m1!=1 && n1!=1)) { - mexErrMsgTxt("COLLINCOMB requires that ROWS be an index vector of type double."); - } - rownum = (m1 > n1) ? m1 : n1; /* the number of rows in the linear combination */ - rownumspecified = 1; - rows = mxGetPr(ROWS_IN); - - m1 = mxGetM(COLS_IN2); - n1 = mxGetN(COLS_IN2); - if (!mxIsDouble(COLS_IN2) || mxIsComplex(COLS_IN2) || (m1!=1 && n1!=1)) { - mexErrMsgTxt("COLLINCOMB requires that COLS be an index vector of type double."); - } - colnum = (m1 > n1) ? m1 : n1; /* the number of columns in the linear combination */ - - m2 = mxGetM(X_IN2); - n2 = mxGetN(X_IN2); - if (!mxIsDouble(X_IN2) || mxIsComplex(X_IN2) || (m2!=1 && n2!=1)) { - mexErrMsgTxt("COLLINCOMB requires that X be a double vector."); - } - - if (m2!=colnum && n2!=colnum) { - mexErrMsgTxt("The length of X does not match the number of columns in COLS."); - } - - Y_OUT = mxCreateDoubleMatrix(rownum, 1, mxREAL); - cols = mxGetPr(COLS_IN2); - x = mxGetPr(X_IN2); - } - - - /* Assign pointers to the various parameters */ - A = mxGetPr(A_IN); - y = mxGetPr(Y_OUT); - - - if (rownumspecified) { - - /* check row indices */ - - row_ids = (mwIndex*)mxMalloc(rownum*sizeof(mwIndex)); - - for (i=0; i<rownum; ++i) { - row_ids[i] = (mwIndex)(rows[i]+0.1)-1; - if (row_ids[i]<0 || row_ids[i]>=m) { - mexErrMsgTxt("Row index in ROWS is out of range."); - } - } - - /* Do the actual computation */ - for (i=0; i<colnum; ++i) { - col_id = (mwIndex)(cols[i]+0.1)-1; - if (col_id<0 || col_id>=n) { - mexErrMsgTxt("Column index in COLS is out of range."); - } - for (j=0; j<rownum; ++j) { - y[j] += A[m*col_id+row_ids[j]]*x[i]; - } - } - - mxFree(row_ids); - } - - else { - - /* Do the actual computation */ - for (i=0; i<colnum; ++i) { - col_id = (mwIndex)(cols[i]+0.1)-1; - if (col_id<0 || col_id>=n) { - mexErrMsgTxt("Column index in COLS is out of range."); - } - for (j=0; j<m; ++j) { - y[j] += A[m*col_id+j]*x[i]; - } - } - } - - return; -}
--- a/Problems/private/collincomb.m Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,27 +0,0 @@ -%COLLINCOMB Linear combination of matrix columns. -% Y = COLLINCOMB(A,COLS,X) computes a linear combination of the columns of -% the matrix A. The column indices are specified in the vector COLS, and -% the correspoinding coefficients are specified in the vector X. The -% vectors COLS and X must be of the same length. -% -% The call Y = COLLINCOMB(A,COLS,X) is essentially equivalent to -% -% Y = A(:,COLS)*X . -% -% However, it is implemented much more efficiently. -% -% Y = COLLINCOMB(A,ROWS,COLS,X) only works on the rows of A specified -% in ROWS, returning a vector of length equal to ROWS. This call is -% essentially equivalent to the command -% -% Y = A(ROWS,COLS)*X . -% -% See also ROWLINCOMB. - - -% Ron Rubinstein -% Computer Science Department -% Technion, Haifa 32000 Israel -% ronrubin@cs -% -% April 2009
--- a/Problems/private/completeOps.m Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,78 +0,0 @@ -function data = completeOps(data) - -% Copyright 2008, Ewout van den Berg and Michael P. Friedlander -% http://www.cs.ubc.ca/labs/scl/sparco -% $Id: completeOps.m 1040 2008-06-26 20:29:02Z ewout78 $ - -operators = {}; -flagM = 0; if isfield(data,'M'), flagM = 1; end; -flagB = 0; if isfield(data,'B'), flagB = 1; end; - -if (~flagM) && (~flagB) - error('At least one of the operators M or B has be to given.'); -end - -% Define measurement matrix if needed -if ~flagM - info = data.B([],0); - data.M = opDirac(info{1}); -else - operators{end+1} = data.M; -end - -% Define sparsity basis if needed -if ~flagB - info = data.M([],0); - data.B = opDirac(info{2}); -else - operators{end+1} = data.B; -end - -% Define operator A if needed -if ~isfield(data,'A') - if (length(operators) > 1) - data.A = opFoG(operators{:}); - else - data.A = operators{1}; - end -end - -% Define empty solution if needed -if ~isfield(data,'x0') - data.x0 = []; -end - -% Define the operator size and string -opInfo = data.A([],0); -data.sizeA = [opInfo{1},opInfo{2}]; -opInfo = data.B([],0); -data.sizeB = [opInfo{1},opInfo{2}]; -opInfo = data.M([],0); -data.sizeM = [opInfo{1},opInfo{2}]; -data.op.strA = opToString(data.A); -data.op.strB = opToString(data.B); -data.op.strM = opToString(data.M); - -% Get the size of the desired signal -if ~isfield(data,'signalSize') - if ~isfield(data,'signal') - error(['At least one of the fields signal ', ... - 'or signalSize must be given.']); - end - data.signalSize = size(data.signal); -end - -% Reconstruct signal from sparse coefficients -if ~isfield(data,'reconstruct') - data.reconstruct = @(x) reshape(data.B(x,1),data.signalSize); -end - -% Reorder the fields (sort alphabetically) -m = fieldnames(data); -m = sort(m); -dataReorder = struct(); -for i=1:length(m) - eval(sprintf('dataReorder.%s = data.%s;',m{i},m{i})); -end - -data = dataReorder;
--- a/Problems/private/countcover.m Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,33 +0,0 @@ -function cnt = countcover(sz,blocksize,stepsize) -%COUNTCOVER Covering of signal samples by blocks -% CNT = COUNTCOVER(SZ,BLOCKSIZE,STEPSIZE) assumes a p-dimensional signal -% of size SZ=[N1 N2 ... Np] covered by (possibly overlapping) blocks of -% size BLOCKSIZE=[M1 M2 ... Mp]. The blocks start at position (1,1,..,1) -% and are shifted between them by steps of size STEPSIZE=[S1 S2 ... Sp]. -% COUNTCOVER returns a matrix the same size as the signal, containing in -% each entry the number of blocks covering that sample. -% -% See also IM2COLSTEP, COL2IMSTEP, IM2COL. - - -% Ron Rubinstein -% Computer Science Department -% Technion, Haifa 32000 Israel -% ronrubin@cs -% -% August 2008 - - -cnt = ones(sz); -for k = 1:length(sz) - - % this code is modified from function NDGRID, so it computes one - % output argument of NDGRID at a time (to conserve memory) - ids = (1:sz(k))'; - s = sz; s(k) = []; - ids = reshape(ids(:,ones(1,prod(s))),[length(ids) s]); - ids = permute(ids,[2:k 1 k+1:length(sz)]); - - cnt = cnt .* max( min(floor((ids-1)/stepsize(k)),floor((sz(k)-blocksize(k))/stepsize(k))) - ... - max(ceil((ids-blocksize(k))/stepsize(k)),0) + 1 , 0 ); -end
--- a/Problems/private/diag_ids.m Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,30 +0,0 @@ -function id = diag_ids(x,k) -%DIAG_IDS Indices of matrix diagonal elements. -% ID = DIAG_IDS(X) returns the indices of the main diagonal of X. -% -% ID = DIAG_IDS(X,K) returns the indices of the K-th diagonal. K=0 -% represents the main diagonal, positive values are above the main -% diagonal and negative values are below the main diagonal. - - -% Ron Rubinstein -% Computer Science Department -% Technion, Haifa 32000 Israel -% ronrubin@cs -% -% September 2006 - - -if (nargin==1), k=0; end - -[m,n] = size(x); -l = max(m,n); - -if (k<=0) - id = (0:l-1)*m + (1:l) - k; -else - id = (0:l-1)*m + (1:l) + k*m; -end - -if (l-k>m), id = id(1:end-(l-k-m)); end -if (l+k>n), id = id(1:end-(l+k-n)); end
--- a/Problems/private/dictdist.m Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,61 +0,0 @@ -function [dist,ratio] = dictdist(approxD,D,epsilon) -%DICTDIST Distance between dictionaries. -% [DIST,RATIO] = DICTDIST(APPROXD,D) computes the distance between the -% approximate dictionary APPROXD and the true dictionary D, where APPROXD -% is NxK and D is NxM. -% -% The distance between the dictionary APPROXD and a single atom A of D is -% defined as: -% -% DIST(APPROXD,A) = min { 1-abs(APPROXD(:,i)' * A) } -% i -% -% The distance between the dictionaries APPROXD and D is defined as: -% -% DIST(APPROXD,D) = sum { dist(APPROXD, D(:,k)) } / M -% k -% -% Note that 0 <= DIST(APPROXD,D) <= 1, where 0 implies that all atoms in D -% appear in APPROXD, and 1 implies that the atoms of D are orthogonal to -% APPROXD. -% -% The similarity ratio between APPROXD and D is defined as: -% -% RATIO(APPROXD,D) = #(atoms in D that appear in APPROXD) / M -% -% where two atoms are considered identical when DIST(A1,A2) < EPSILON with -% EPSILON=0.01 by default. Note that 0 <= RATIO(APPROXD,D) <= 1, where 0 -% means APPROXD and D have no identical atoms, and 1 means that all atoms -% of D appear in APPROXD. -% -% [DIST,RATIO] = DICTDIST(DICT1,DICT2,EPSILON) specifies a different value -% for EPSILON. - -% Ron Rubinstein -% Computer Science Department -% Technion, Haifa 32000 Israel -% ronrubin@cs -% -% October 2007 - - -if (nargin < 3), epsilon = 0.01; end - -[n,m] = size(D); - -approxD = normcols(approxD*spdiag(sign(approxD(1,:)))); -D = normcols(D*spdiag(sign(D(1,:)))); - -identical_atoms = 0; -dist = 0; - -for i = 1:m - atom = D(:,i); - distances = 1-abs(atom'*approxD); - mindist = min(distances); - dist = dist + mindist; - identical_atoms = identical_atoms + (mindist < epsilon); -end - -dist = dist / m; -ratio = identical_atoms / m;
--- a/Problems/private/genPDF.m Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,113 +0,0 @@ -function [pdf,val] = genPDF(imSize,p,pctg,distType,radius,disp) - -%[pdf,val] = genPDF(imSize,p,pctg [,distType,radius,disp]) -% -% generates a pdf for a 1d or 2d random sampling pattern -% with polynomial variable density sampling -% -% Input: -% imSize - size of matrix or vector -% p - power of polynomial -% pctg - partial sampling factor e.g. 0.5 for half -% distType - 1 or 2 for L1 or L2 distance measure -% radius - radius of fully sampled center -% disp - display output -% -% Output: -% pdf - the pdf -% val - min sampling density -% -% -% Example: -% [pdf,val] = genPDF([128,128],2,0.5,2,0,1); -% -% (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: genPDF.m 1040 2008-06-26 20:29:02Z ewout78 $ - - -if nargin < 4 - distType = 2; -end - -if nargin < 5 - radius = 0; -end - -if nargin < 6 - disp = 0; -end - - -minval=0; -maxval=1; -val = 0.5; - -if length(imSize)==1 - imSize = [imSize,1]; -end - -sx = imSize(1); -sy = imSize(2); -PCTG = floor(pctg*sx*sy); - - -if sum(imSize==1)==0 % 2D - [x,y] = meshgrid(linspace(-1,1,sy),linspace(-1,1,sx)); - switch distType - case 1 - r = max(abs(x),abs(y)); - otherwise - r = sqrt(x.^2+y.^2); - r = r/max(abs(r(:))); - end - -else %1d - r = abs(linspace(-1,1,max(sx,sy))); -end - - - - -idx = find(r<radius); - -pdf = (1-r).^p; pdf(idx) = 1; -if floor(sum(pdf(:))) > PCTG - error('infeasible without undersampling dc, increase p'); -end - -% begin bisection -while(1) - val = minval/2 + maxval/2; - pdf = (1-r).^p + val; pdf(find(pdf>1)) = 1; pdf(idx)=1; - N = floor(sum(pdf(:))); - if N > PCTG % infeasible - maxval=val; - end - if N < PCTG % feasible, but not optimal - minval=val; - end - if N==PCTG % optimal - break; - end -end - -if disp - figure, - subplot(211), imshow(pdf) - if sum(imSize==1)==0 - subplot(212), plot(pdf(end/2+1,:)); - else - subplot(212), plot(pdf); - end -end - - - - - -
--- a/Problems/private/genSampling.m Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,53 +0,0 @@ -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))<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); - -
--- a/Problems/private/im2colstep.c Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,131 +0,0 @@ -/************************************************************************** - * - * File name: im2colstep.c - * - * Ron Rubinstein - * Computer Science Department - * Technion, Haifa 32000 Israel - * ronrubin@cs - * - * Last Updated: 31.8.2009 - * - *************************************************************************/ - - -#include "mex.h" -#include <string.h> - - -/* Input Arguments */ - -#define X_IN prhs[0] -#define SZ_IN prhs[1] -#define S_IN prhs[2] - - -/* Output Arguments */ - -#define B_OUT plhs[0] - - -void mexFunction(int nlhs, mxArray *plhs[], - int nrhs, const mxArray*prhs[]) - -{ - double *x, *b, *s; - mwSize sz[3], stepsize[3], n[3], ndims; - mwIndex i, j, k, l, m, blocknum; - - - /* Check for proper number of arguments */ - - if (nrhs < 2 || nrhs > 3) { - mexErrMsgTxt("Invalid number of input arguments."); - } else if (nlhs > 1) { - mexErrMsgTxt("Too many output arguments."); - } - - - /* Check the the input dimensions */ - - ndims = mxGetNumberOfDimensions(X_IN); - - if (!mxIsDouble(X_IN) || mxIsComplex(X_IN) || ndims>3) { - mexErrMsgTxt("X should be a 2-D or 3-D double matrix."); - } - if (!mxIsDouble(SZ_IN) || mxIsComplex(SZ_IN) || mxGetNumberOfDimensions(SZ_IN)>2 || mxGetM(SZ_IN)*mxGetN(SZ_IN)!=ndims) { - mexErrMsgTxt("Invalid block size."); - } - if (nrhs == 3) { - if (!mxIsDouble(S_IN) || mxIsComplex(S_IN) || mxGetNumberOfDimensions(S_IN)>2 || mxGetM(S_IN)*mxGetN(S_IN)!=ndims) { - mexErrMsgTxt("Invalid step size."); - } - } - - - /* Get parameters */ - - s = mxGetPr(SZ_IN); - if (s[0]<1 || s[1]<1 || (ndims==3 && s[2]<1)) { - mexErrMsgTxt("Invalid block size."); - } - sz[0] = (mwSize)(s[0] + 0.01); - sz[1] = (mwSize)(s[1] + 0.01); - sz[2] = ndims==3 ? (mwSize)(s[2] + 0.01) : 1; - - if (nrhs == 3) { - s = mxGetPr(S_IN); - if (s[0]<1 || s[1]<1 || (ndims==3 && s[2]<1)) { - mexErrMsgTxt("Invalid step size."); - } - stepsize[0] = (mwSize)(s[0] + 0.01); - stepsize[1] = (mwSize)(s[1] + 0.01); - stepsize[2] = ndims==3 ? (mwSize)(s[2] + 0.01) : 1; - } - else { - stepsize[0] = stepsize[1] = stepsize[2] = 1; - } - - n[0] = (mxGetDimensions(X_IN))[0]; - n[1] = (mxGetDimensions(X_IN))[1]; - n[2] = ndims==3 ? (mxGetDimensions(X_IN))[2] : 1; - - if (n[0]<sz[0] || n[1]<sz[1] || (ndims==3 && n[2]<sz[2])) { - mexErrMsgTxt("Block size too large."); - } - - - /* Create a matrix for the return argument */ - - B_OUT = mxCreateDoubleMatrix(sz[0]*sz[1]*sz[2], ((n[0]-sz[0])/stepsize[0]+1)*((n[1]-sz[1])/stepsize[1]+1)*((n[2]-sz[2])/stepsize[2]+1), mxREAL); - - - /* Assign pointers */ - - x = mxGetPr(X_IN); - b = mxGetPr(B_OUT); - - - /* Do the actual computation */ - - blocknum = 0; - - /* iterate over all blocks */ - for (k=0; k<=n[2]-sz[2]; k+=stepsize[2]) { - for (j=0; j<=n[1]-sz[1]; j+=stepsize[1]) { - for (i=0; i<=n[0]-sz[0]; i+=stepsize[0]) { - - /* copy single block */ - for (m=0; m<sz[2]; m++) { - for (l=0; l<sz[1]; l++) { - memcpy(b + blocknum*sz[0]*sz[1]*sz[2] + m*sz[0]*sz[1] + l*sz[0], x+(k+m)*n[0]*n[1]+(j+l)*n[0]+i, sz[0]*sizeof(double)); - } - } - blocknum++; - - } - } - } - - return; -}
--- a/Problems/private/im2colstep.m Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,31 +0,0 @@ -%IM2COLSTEP Rearrange matrix blocks into columns. -% B = IM2COLSTEP(A,[N1 N2]) converts each sliding N1-by-N2 block of the -% 2-D matrix A into a column of B, with no zero padding. B has N1*N2 rows -% and will contain as many columns as there are N1-by-N2 neighborhoods in -% A. Each column of B contains a neighborhood of A reshaped as NHOOD(:), -% where NHOOD is a matrix containing an N1-by-N2 neighborhood of A. -% -% B = IM2COLSTEP(A,[N1 N2],[S1 S2]) extracts neighborhoods of A with a -% step size of (S1,S2) between them. The first extracted neighborhood is -% A(1:N1,1:N2), and the rest are of the form A((1:N1)+i*S1,(1:N2)+j*S2). -% Note that to ensure coverage of all A by neighborhoods, -% (size(A,i)-Ni)/Si must be whole for i=1,2. The default function behavior -% corresponds to [S1 S2] = [1 1]. Setting S1>=N1 and S2>=N2 results in no -% overlap between the neighborhoods. -% -% B = IM2COLSTEP(A,[N1 N2 N3],[S1 S2 S3]) operates on a 3-D matrix A. The -% step size [S1 S2 S3] may be ommitted, and defaults to [1 1 1]. -% -% Note: the call IM2COLSTEP(A,[N1 N2]) produces the same output as -% Matlab's IM2COL(A,[N1 N2],'sliding'). However, it is significantly -% faster. -% -% See also COL2IMSTEP, IM2COL, COUNTCOVER. - - -% Ron Rubinstein -% Computer Science Department -% Technion, Haifa 32000 Israel -% ronrubin@cs -% -% August 2009
--- a/Problems/private/imnormalize.m Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,19 +0,0 @@ -function y = imnormalize(x) -%IMNORMALIZE Normalize image values. -% Y = IMNORMALIZE(X) linearly transforms the intensity values of the image -% X to tightly cover the range [0,1]. If X has more than one channel, the -% channels are handled as one and normalized using the same transform. - - -% Ron Rubinstein -% Computer Science Department -% Technion, Haifa 32000 Israel -% ronrubin@cs -% -% May 2004 - - -maxval = max(x(:)); -minval = min(x(:)); - -y = (x-minval) / (maxval-minval);
--- a/Problems/private/iswhole.m Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,17 +0,0 @@ -function z = iswhole(x,epsilon) -%ISWHOLE Determine whole numbers with finite precision. -% Z = ISWHOLE(X,EPSILON) returns a matrix the same size as X, containing -% 1's where X is whole up to precision EPSILON, and 0's elsewhere. -% -% Z = ISWHOLE(X) uses the default value of EPSILON=1e-6. - - -% Ron Rubinstein -% Computer Science Department -% Technion, Haifa 32000 Israel -% ronrubin@cs -% -% May 2005 - -if (nargin<2), epsilon = 1e-6; end -z = abs(round(x)-x) < epsilon;
--- a/Problems/private/make.m Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,40 +0,0 @@ -function make -%MAKE Build the KSVDBox MEX support files. -% MAKE compiles the KSVDBox supporting MEX functions, using Matlab's -% default MEX compiler. If the MEX compiler has not been set-up before, -% please run -% -% mex -setup -% -% before using this MAKE file. - -% Ron Rubinstein -% Computer Science Department -% Technion, Haifa 32000 Israel -% ronrubin@cs -% -% August 2009 - - -% detect platform - -compstr = computer; -is64bit = strcmp(compstr(end-1:end),'64'); - - -% compilation parameters - -compile_params = cell(0); -if (is64bit) - compile_params{1} = '-largeArrayDims'; -end - - -% Compile files % - -sourcefiles = {{'addtocols.c'}, {'collincomb.c'}, {'rowlincomb.c'}, {'sprow.c','mexutils.c'}, {'im2colstep.c'}, {'col2imstep.c'}}; - -for i = 1:length(sourcefiles) - printf('Compiling %s...', sourcefiles{i}{1}); - mex(sourcefiles{i}{:},compile_params{:}); -end
--- a/Problems/private/mexutils.c Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,79 +0,0 @@ -/************************************************************************** - * - * File name: mexutils.c - * - * Ron Rubinstein - * Computer Science Department - * Technion, Haifa 32000 Israel - * ronrubin@cs - * - * Last Updated: 15.8.2009 - * - *************************************************************************/ - -#include "mexutils.h" -#include <math.h> - - - -/* verify that the mxArray contains a double matrix */ - -void checkmatrix(const mxArray *param, char *fname, char *pname) -{ - char errmsg[100]; - sprintf(errmsg, "%.15s requires that %.25s be a double matrix.", fname, pname); - if (!mxIsDouble(param) || mxIsComplex(param) || mxGetNumberOfDimensions(param)>2) { - mexErrMsgTxt(errmsg); - } -} - - -/* verify that the mxArray contains a 1-D double vector */ - -void checkvector(const mxArray *param, char *fname, char *pname) -{ - char errmsg[100]; - sprintf(errmsg, "%.15s requires that %.25s be a double vector.", fname, pname); - if (!mxIsDouble(param) || mxIsComplex(param) || mxGetNumberOfDimensions(param)>2 || (mxGetM(param)!=1 && mxGetN(param)!=1)) { - mexErrMsgTxt(errmsg); - } -} - - -/* verify that the mxArray contains a double scalar */ - -void checkscalar(const mxArray *param, char *fname, char *pname) -{ - char errmsg[100]; - sprintf(errmsg, "%.15s requires that %.25s be a double scalar.", fname, pname); - if (!mxIsDouble(param) || mxIsComplex(param) || mxGetNumberOfDimensions(param)>2 || - mxGetM(param)!=1 || mxGetN(param)!=1) - { - mexErrMsgTxt(errmsg); - } -} - - -/* verify that the mxArray contains a sparse matrix */ - -void checksparse(const mxArray *param, char *fname, char *pname) -{ - char errmsg[100]; - sprintf(errmsg, "%.15s requires that %.25s be sparse.", fname, pname); - if (!mxIsSparse(param)) { - mexErrMsgTxt(errmsg); - } -} - - -/* verify that the mxArray contains a 1-dimensional cell array */ - -void checkcell_1d(const mxArray *param, char *fname, char *pname) -{ - char errmsg[100]; - sprintf(errmsg, "%.15s requires that %.25s be a 1-D cell array.", fname, pname); - if (!mxIsCell(param) || mxGetNumberOfDimensions(param)>2 || (mxGetM(param)!=1 && mxGetN(param)!=1)) { - mexErrMsgTxt(errmsg); - } -} -
--- a/Problems/private/mexutils.h Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,103 +0,0 @@ -/************************************************************************** - * - * File name: mexutils.h - * - * Ron Rubinstein - * Computer Science Department - * Technion, Haifa 32000 Israel - * ronrubin@cs - * - * Last Updated: 18.8.2009 - * - * Utility functions for MEX files. - * - *************************************************************************/ - - -#ifndef __MEX_UTILS_H__ -#define __MEX_UTILS_H__ - -#include "mex.h" - - - -/************************************************************************** - * Function checkmatrix: - * - * Verify that the specified mxArray is real, of type double, and has - * no more than two dimensions. If not, an error message is printed - * and the mex file terminates. - * - * Parameters: - * param - the mxArray to be checked - * fname - the name of the function where the error occured (15 characters or less) - * pname - the name of the parameter (25 characters or less) - * - **************************************************************************/ -void checkmatrix(const mxArray *param, char *fname, char *pname); - - -/************************************************************************** - * Function checkvector: - * - * Verify that the specified mxArray is 1-D, real, and of type double. The - * vector may be a column or row vector. Otherwise, an error message is - * printed and the mex file terminates. - * - * Parameters: - * param - the mxArray to be checked - * fname - the name of the function where the error occured (15 characters or less) - * pname - the name of the parameter (25 characters or less) - * - **************************************************************************/ -void checkvector(const mxArray *param, char *fname, char *pname); - - -/************************************************************************** - * Function checkscalar: - * - * Verify that the specified mxArray represents a real double scalar value. - * If not, an error message is printed and the mex file terminates. - * - * Parameters: - * param - the mxArray to be checked - * fname - the name of the function where the error occured (15 characters or less) - * pname - the name of the parameter (25 characters or less) - * - **************************************************************************/ -void checkscalar(const mxArray *param, char *fname, char *pname); - - -/************************************************************************** - * Function checksparse: - * - * Verify that the specified mxArray contains a sparse matrix. If not, - * an error message is printed and the mex file terminates. - * - * Parameters: - * param - the mxArray to be checked - * fname - the name of the function where the error occured (15 characters or less) - * pname - the name of the parameter (25 characters or less) - * - **************************************************************************/ -void checksparse(const mxArray *param, char *fname, char *pname); - - -/************************************************************************** - * Function checkcell_1d: - * - * Verify that the specified mxArray is a 1-D cell array. The cell array - * may be arranged as either a column or a row. If not, an error message - * is printed and the mex file terminates. - * - * Parameters: - * param - the mxArray to be checked - * fname - the name of the function where the error occured (15 characters or less) - * pname - the name of the parameter (25 characters or less) - * - **************************************************************************/ -void checkcell_1d(const mxArray *param, char *fname, char *pname); - - -#endif -
--- a/Problems/private/myblas.c Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,663 +0,0 @@ -/************************************************************************** - * - * 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); -} -*/ - - - -
--- a/Problems/private/myblas.h Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,492 +0,0 @@ -/************************************************************************** - * - * 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 -
--- a/Problems/private/normcols.m Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,17 +0,0 @@ -function y = normcols(x) -%NORMCOLS Normalize matrix columns. -% Y = NORMCOLS(X) normalizes the columns of X to unit length, returning -% the result as Y. -% -% See also ADDTOCOLS. - - -% Ron Rubinstein -% Computer Science Department -% Technion, Haifa 32000 Israel -% ronrubin@cs -% -% April 2009 - - -y = x*spdiag(1./sqrt(sum(x.*x)));
--- a/Problems/private/omp2mex.c Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,156 +0,0 @@ -/************************************************************************** - * - * 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; -}
--- a/Problems/private/omp2mex.m Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,23 +0,0 @@ -%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
--- a/Problems/private/ompcore.c Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,409 +0,0 @@ -/************************************************************************** - * - * File name: ompcore.c - * - * Ron Rubinstein - * Computer Science Department - * Technion, Haifa 32000 Israel - * ronrubin@cs - * - * Last Updated: 25.8.2009 - * - *************************************************************************/ - - -#include "ompcore.h" -#include "omputils.h" -#include "ompprof.h" -#include "myblas.h" -#include <math.h> -#include <string.h> - - - -/****************************************************************************** - * * - * Batch-OMP Implementation * - * * - ******************************************************************************/ - -mxArray* ompcore(double D[], double x[], double DtX[], double XtX[], double G[], mwSize n, mwSize m, mwSize L, - int T, double eps, int gamma_mode, int profile, double msg_delta, int erroromp) -{ - - profdata pd; - mxArray *Gamma; - mwIndex i, j, signum, pos, *ind, *gammaIr, *gammaJc, gamma_count; - mwSize allocated_coefs, allocated_cols; - int DtX_specified, XtX_specified, batchomp, standardomp, *selected_atoms; - double *alpha, *r, *Lchol, *c, *Gsub, *Dsub, sum, *gammaPr, *tempvec1, *tempvec2; - double eps2, resnorm, delta, deltaprev, secs_remain; - int mins_remain, hrs_remain; - clock_t lastprint_time, starttime; - - - - /*** status flags ***/ - - DtX_specified = (DtX!=0); /* indicates whether D'*x was provided */ - XtX_specified = (XtX!=0); /* indicates whether sum(x.*x) was provided */ - - standardomp = (G==0); /* batch-omp or standard omp are selected depending on availability of G */ - batchomp = !standardomp; - - - - /*** allocate output matrix ***/ - - - if (gamma_mode == FULL_GAMMA) { - - /* allocate full matrix of size m X L */ - - Gamma = mxCreateDoubleMatrix(m, L, mxREAL); - gammaPr = mxGetPr(Gamma); - gammaIr = 0; - gammaJc = 0; - } - else { - - /* allocate sparse matrix with room for allocated_coefs nonzeros */ - - /* for error-omp, begin with L*sqrt(n)/2 allocated nonzeros, otherwise allocate L*T nonzeros */ - allocated_coefs = erroromp ? (mwSize)(ceil(L*sqrt((double)n)/2.0) + 1.01) : L*T; - Gamma = mxCreateSparse(m, L, allocated_coefs, mxREAL); - gammaPr = mxGetPr(Gamma); - gammaIr = mxGetIr(Gamma); - gammaJc = mxGetJc(Gamma); - gamma_count = 0; - gammaJc[0] = 0; - } - - - /*** helper arrays ***/ - - alpha = (double*)mxMalloc(m*sizeof(double)); /* contains D'*residual */ - ind = (mwIndex*)mxMalloc(n*sizeof(mwIndex)); /* indices of selected atoms */ - selected_atoms = (int*)mxMalloc(m*sizeof(int)); /* binary array with 1's for selected atoms */ - c = (double*)mxMalloc(n*sizeof(double)); /* orthogonal projection result */ - - /* current number of columns in Dsub / Gsub / Lchol */ - allocated_cols = erroromp ? (mwSize)(ceil(sqrt((double)n)/2.0) + 1.01) : T; - - /* Cholesky decomposition of D_I'*D_I */ - Lchol = (double*)mxMalloc(n*allocated_cols*sizeof(double)); - - /* temporary vectors for various computations */ - tempvec1 = (double*)mxMalloc(m*sizeof(double)); - tempvec2 = (double*)mxMalloc(m*sizeof(double)); - - if (batchomp) { - /* matrix containing G(:,ind) - the columns of G corresponding to the selected atoms, in order of selection */ - Gsub = (double*)mxMalloc(m*allocated_cols*sizeof(double)); - } - else { - /* matrix containing D(:,ind) - the selected atoms from D, in order of selection */ - Dsub = (double*)mxMalloc(n*allocated_cols*sizeof(double)); - - /* stores the residual */ - r = (double*)mxMalloc(n*sizeof(double)); - } - - if (!DtX_specified) { - /* contains D'*x for the current signal */ - DtX = (double*)mxMalloc(m*sizeof(double)); - } - - - - /*** initializations for error omp ***/ - - if (erroromp) { - eps2 = eps*eps; /* compute eps^2 */ - if (T<0 || T>n) { /* unspecified max atom num - set max atoms to n */ - T = n; - } - } - - - - /*** initialize timers ***/ - - initprofdata(&pd); /* initialize profiling counters */ - starttime = clock(); /* record starting time for eta computations */ - lastprint_time = starttime; /* time of last status display */ - - - - /********************** perform omp for each signal **********************/ - - - - for (signum=0; signum<L; ++signum) { - - - /* initialize residual norm and deltaprev for error-omp */ - - if (erroromp) { - if (XtX_specified) { - resnorm = XtX[signum]; - } - else { - resnorm = dotprod(x+n*signum, x+n*signum, n); - addproftime(&pd, XtX_TIME); - } - deltaprev = 0; /* delta tracks the value of gamma'*G*gamma */ - } - else { - /* ignore residual norm stopping criterion */ - eps2 = 0; - resnorm = 1; - } - - - if (resnorm>eps2 && T>0) { - - /* compute DtX */ - - if (!DtX_specified) { - matT_vec(1, D, x+n*signum, DtX, n, m); - addproftime(&pd, DtX_TIME); - } - - - /* initialize alpha := DtX */ - - memcpy(alpha, DtX + m*signum*DtX_specified, m*sizeof(double)); - - - /* mark all atoms as unselected */ - - for (i=0; i<m; ++i) { - selected_atoms[i] = 0; - } - - } - - - /* main loop */ - - i=0; - while (resnorm>eps2 && i<T) { - - /* index of next atom */ - - pos = maxabs(alpha, m); - addproftime(&pd, MAXABS_TIME); - - - /* stop criterion: selected same atom twice, or inner product too small */ - - if (selected_atoms[pos] || alpha[pos]*alpha[pos]<1e-14) { - break; - } - - - /* mark selected atom */ - - ind[i] = pos; - selected_atoms[pos] = 1; - - - /* matrix reallocation */ - - if (erroromp && i>=allocated_cols) { - - allocated_cols = (mwSize)(ceil(allocated_cols*MAT_INC_FACTOR) + 1.01); - - Lchol = (double*)mxRealloc(Lchol,n*allocated_cols*sizeof(double)); - - batchomp ? (Gsub = (double*)mxRealloc(Gsub,m*allocated_cols*sizeof(double))) : - (Dsub = (double*)mxRealloc(Dsub,n*allocated_cols*sizeof(double))) ; - } - - - /* append column to Gsub or Dsub */ - - if (batchomp) { - memcpy(Gsub+i*m, G+pos*m, m*sizeof(double)); - } - else { - memcpy(Dsub+i*n, D+pos*n, n*sizeof(double)); - } - - - /*** Cholesky update ***/ - - if (i==0) { - *Lchol = 1; - } - else { - - /* incremental Cholesky decomposition: compute next row of Lchol */ - - if (standardomp) { - matT_vec(1, Dsub, D+n*pos, tempvec1, n, i); /* compute tempvec1 := Dsub'*d where d is new atom */ - addproftime(&pd, DtD_TIME); - } - else { - vec_assign(tempvec1, Gsub+i*m, ind, i); /* extract tempvec1 := Gsub(ind,i) */ - } - backsubst('L', Lchol, tempvec1, tempvec2, n, i); /* compute tempvec2 = Lchol \ tempvec1 */ - for (j=0; j<i; ++j) { /* write tempvec2 to end of Lchol */ - Lchol[j*n+i] = tempvec2[j]; - } - - /* compute Lchol(i,i) */ - sum = 0; - for (j=0; j<i; ++j) { /* compute sum of squares of last row without Lchol(i,i) */ - sum += SQR(Lchol[j*n+i]); - } - if ( (1-sum) <= 1e-14 ) { /* Lchol(i,i) is zero => selected atoms are dependent */ - break; - } - Lchol[i*n+i] = sqrt(1-sum); - } - - addproftime(&pd, LCHOL_TIME); - - i++; - - - /* perform orthogonal projection and compute sparse coefficients */ - - vec_assign(tempvec1, DtX + m*signum*DtX_specified, ind, i); /* extract tempvec1 = DtX(ind) */ - cholsolve('L', Lchol, tempvec1, c, n, i); /* solve LL'c = tempvec1 for c */ - addproftime(&pd, COMPCOEF_TIME); - - - /* update alpha = D'*residual */ - - if (standardomp) { - mat_vec(-1, Dsub, c, r, n, i); /* compute r := -Dsub*c */ - vec_sum(1, x+n*signum, r, n); /* compute r := x+r */ - - - /*memcpy(r, x+n*signum, n*sizeof(double)); /* assign r := x */ - /*mat_vec1(-1, Dsub, c, 1, r, n, i); /* compute r := r-Dsub*c */ - - addproftime(&pd, COMPRES_TIME); - matT_vec(1, D, r, alpha, n, m); /* compute alpha := D'*r */ - addproftime(&pd, DtR_TIME); - - /* update residual norm */ - if (erroromp) { - resnorm = dotprod(r, r, n); - addproftime(&pd, UPDATE_RESNORM_TIME); - } - } - else { - mat_vec(1, Gsub, c, tempvec1, m, i); /* compute tempvec1 := Gsub*c */ - memcpy(alpha, DtX + m*signum*DtX_specified, m*sizeof(double)); /* set alpha = D'*x */ - vec_sum(-1, tempvec1, alpha, m); /* compute alpha := alpha - tempvec1 */ - addproftime(&pd, UPDATE_DtR_TIME); - - /* update residual norm */ - if (erroromp) { - vec_assign(tempvec2, tempvec1, ind, i); /* assign tempvec2 := tempvec1(ind) */ - delta = dotprod(c,tempvec2,i); /* compute c'*tempvec2 */ - resnorm = resnorm - delta + deltaprev; /* residual norm update */ - deltaprev = delta; - addproftime(&pd, UPDATE_RESNORM_TIME); - } - } - } - - - /*** generate output vector gamma ***/ - - if (gamma_mode == FULL_GAMMA) { /* write the coefs in c to their correct positions in gamma */ - for (j=0; j<i; ++j) { - gammaPr[m*signum + ind[j]] = c[j]; - } - } - else { - /* sort the coefs by index before writing them to gamma */ - quicksort(ind,c,i); - addproftime(&pd, INDEXSORT_TIME); - - /* gamma is full - reallocate */ - if (gamma_count+i >= allocated_coefs) { - - while(gamma_count+i >= allocated_coefs) { - allocated_coefs = (mwSize)(ceil(GAMMA_INC_FACTOR*allocated_coefs) + 1.01); - } - - mxSetNzmax(Gamma, allocated_coefs); - mxSetPr(Gamma, mxRealloc(gammaPr, allocated_coefs*sizeof(double))); - mxSetIr(Gamma, mxRealloc(gammaIr, allocated_coefs*sizeof(mwIndex))); - - gammaPr = mxGetPr(Gamma); - gammaIr = mxGetIr(Gamma); - } - - /* append coefs to gamma and update the indices */ - for (j=0; j<i; ++j) { - gammaPr[gamma_count] = c[j]; - gammaIr[gamma_count] = ind[j]; - gamma_count++; - } - gammaJc[signum+1] = gammaJc[signum] + i; - } - - - - /*** display status messages ***/ - - if (msg_delta>0 && (clock()-lastprint_time)/(double)CLOCKS_PER_SEC >= msg_delta) - { - lastprint_time = clock(); - - /* estimated remainig time */ - secs2hms( ((L-signum-1)/(double)(signum+1)) * ((lastprint_time-starttime)/(double)CLOCKS_PER_SEC) , - &hrs_remain, &mins_remain, &secs_remain); - - mexPrintf("omp: signal %d / %d, estimated remaining time: %02d:%02d:%05.2f\n", - signum+1, L, hrs_remain, mins_remain, secs_remain); - mexEvalString("drawnow;"); - } - - } - - /* end omp */ - - - - /*** print final messages ***/ - - if (msg_delta>0) { - mexPrintf("omp: signal %d / %d\n", signum, L); - } - - if (profile) { - printprofinfo(&pd, erroromp, batchomp, L); - } - - - - /* free memory */ - - if (!DtX_specified) { - mxFree(DtX); - } - if (standardomp) { - mxFree(r); - mxFree(Dsub); - } - else { - mxFree(Gsub); - } - mxFree(tempvec2); - mxFree(tempvec1); - mxFree(Lchol); - mxFree(c); - mxFree(selected_atoms); - mxFree(ind); - mxFree(alpha); - - return Gamma; -}
--- a/Problems/private/ompcore.h Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,80 +0,0 @@ -/************************************************************************** - * - * 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
--- a/Problems/private/ompmex.c Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,133 +0,0 @@ -/************************************************************************** - * - * 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; -} -
--- a/Problems/private/ompmex.m Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,22 +0,0 @@ -%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
--- a/Problems/private/ompprof.c Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,113 +0,0 @@ -/************************************************************************** - * - * 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); -} -
--- a/Problems/private/ompprof.h Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,106 +0,0 @@ -/************************************************************************** - * - * 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 -
--- a/Problems/private/omputils.c Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,89 +0,0 @@ -/************************************************************************** - * - * 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--; - } - } -}
--- a/Problems/private/omputils.h Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,77 +0,0 @@ -/************************************************************************** - * - * 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 -
--- a/Problems/private/printf.m Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,26 +0,0 @@ -function str = printf(varargin) -%PRINTF Print formatted text to screen. -% PRINTF(FMT,VAL1,VAL2,...) formats the data in VAL1,VAL2,... according to -% the format string FMT, and prints the result to the screen. -% -% The call to PRINTF(FMT,VAL1,VAL2,...) simply invokes the call -% DISP(SPRINTF(FMT,VAL1,VAL2,...)). For a complete description of the -% format string options see function SPRINTF. -% -% STR = PRINTF(...) also returns the formatted string. - - -% Ron Rubinstein -% Computer Science Department -% Technion, Haifa 32000 Israel -% ronrubin@cs -% -% April 2008 - - -if (nargout>0) - str = sprintf(varargin{:}); - disp(str); -else - disp(sprintf(varargin{:})); -end
--- a/Problems/private/private/make.m Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,41 +0,0 @@ -function make -%MAKE Build the OMPBox package. -% MAKE compiles all OMPBox MEX functions, using Matlab's default MEX -% compiler. If the MEX compiler has not been set-up before, please run -% -% mex -setup -% -% before using this MAKE file. - -% Ron Rubinstein -% Computer Science Department -% Technion, Haifa 32000 Israel -% ronrubin@cs -% -% August 2009 - - -% detect platform - -compstr = computer; -is64bit = strcmp(compstr(end-1:end),'64'); - - -% compilation parameters - -compile_params = cell(0); -if (is64bit) - compile_params{1} = '-largeArrayDims'; -end - - -% Compile files % - -ompsources = {'mexutils.c','ompcore.c','omputils.c','myblas.c','ompprof.c'}; - -disp('Compiling ompmex...'); -mex('ompmex.c', ompsources{:},compile_params{:}); - -disp('Compiling omp2mex...'); -mex('omp2mex.c',ompsources{:},compile_params{:}); -
--- a/Problems/private/private/mexutils.c Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,79 +0,0 @@ -/************************************************************************** - * - * File name: mexutils.c - * - * Ron Rubinstein - * Computer Science Department - * Technion, Haifa 32000 Israel - * ronrubin@cs - * - * Last Updated: 15.8.2009 - * - *************************************************************************/ - -#include "mexutils.h" -#include <math.h> - - - -/* verify that the mxArray contains a double matrix */ - -void checkmatrix(const mxArray *param, char *fname, char *pname) -{ - char errmsg[100]; - sprintf(errmsg, "%.15s requires that %.25s be a double matrix.", fname, pname); - if (!mxIsDouble(param) || mxIsComplex(param) || mxGetNumberOfDimensions(param)>2) { - mexErrMsgTxt(errmsg); - } -} - - -/* verify that the mxArray contains a 1-D double vector */ - -void checkvector(const mxArray *param, char *fname, char *pname) -{ - char errmsg[100]; - sprintf(errmsg, "%.15s requires that %.25s be a double vector.", fname, pname); - if (!mxIsDouble(param) || mxIsComplex(param) || mxGetNumberOfDimensions(param)>2 || (mxGetM(param)!=1 && mxGetN(param)!=1)) { - mexErrMsgTxt(errmsg); - } -} - - -/* verify that the mxArray contains a double scalar */ - -void checkscalar(const mxArray *param, char *fname, char *pname) -{ - char errmsg[100]; - sprintf(errmsg, "%.15s requires that %.25s be a double scalar.", fname, pname); - if (!mxIsDouble(param) || mxIsComplex(param) || mxGetNumberOfDimensions(param)>2 || - mxGetM(param)!=1 || mxGetN(param)!=1) - { - mexErrMsgTxt(errmsg); - } -} - - -/* verify that the mxArray contains a sparse matrix */ - -void checksparse(const mxArray *param, char *fname, char *pname) -{ - char errmsg[100]; - sprintf(errmsg, "%.15s requires that %.25s be sparse.", fname, pname); - if (!mxIsSparse(param)) { - mexErrMsgTxt(errmsg); - } -} - - -/* verify that the mxArray contains a 1-dimensional cell array */ - -void checkcell_1d(const mxArray *param, char *fname, char *pname) -{ - char errmsg[100]; - sprintf(errmsg, "%.15s requires that %.25s be a 1-D cell array.", fname, pname); - if (!mxIsCell(param) || mxGetNumberOfDimensions(param)>2 || (mxGetM(param)!=1 && mxGetN(param)!=1)) { - mexErrMsgTxt(errmsg); - } -} -
--- a/Problems/private/private/mexutils.h Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,103 +0,0 @@ -/************************************************************************** - * - * File name: mexutils.h - * - * Ron Rubinstein - * Computer Science Department - * Technion, Haifa 32000 Israel - * ronrubin@cs - * - * Last Updated: 18.8.2009 - * - * Utility functions for MEX files. - * - *************************************************************************/ - - -#ifndef __MEX_UTILS_H__ -#define __MEX_UTILS_H__ - -#include "mex.h" - - - -/************************************************************************** - * Function checkmatrix: - * - * Verify that the specified mxArray is real, of type double, and has - * no more than two dimensions. If not, an error message is printed - * and the mex file terminates. - * - * Parameters: - * param - the mxArray to be checked - * fname - the name of the function where the error occured (15 characters or less) - * pname - the name of the parameter (25 characters or less) - * - **************************************************************************/ -void checkmatrix(const mxArray *param, char *fname, char *pname); - - -/************************************************************************** - * Function checkvector: - * - * Verify that the specified mxArray is 1-D, real, and of type double. The - * vector may be a column or row vector. Otherwise, an error message is - * printed and the mex file terminates. - * - * Parameters: - * param - the mxArray to be checked - * fname - the name of the function where the error occured (15 characters or less) - * pname - the name of the parameter (25 characters or less) - * - **************************************************************************/ -void checkvector(const mxArray *param, char *fname, char *pname); - - -/************************************************************************** - * Function checkscalar: - * - * Verify that the specified mxArray represents a real double scalar value. - * If not, an error message is printed and the mex file terminates. - * - * Parameters: - * param - the mxArray to be checked - * fname - the name of the function where the error occured (15 characters or less) - * pname - the name of the parameter (25 characters or less) - * - **************************************************************************/ -void checkscalar(const mxArray *param, char *fname, char *pname); - - -/************************************************************************** - * Function checksparse: - * - * Verify that the specified mxArray contains a sparse matrix. If not, - * an error message is printed and the mex file terminates. - * - * Parameters: - * param - the mxArray to be checked - * fname - the name of the function where the error occured (15 characters or less) - * pname - the name of the parameter (25 characters or less) - * - **************************************************************************/ -void checksparse(const mxArray *param, char *fname, char *pname); - - -/************************************************************************** - * Function checkcell_1d: - * - * Verify that the specified mxArray is a 1-D cell array. The cell array - * may be arranged as either a column or a row. If not, an error message - * is printed and the mex file terminates. - * - * Parameters: - * param - the mxArray to be checked - * fname - the name of the function where the error occured (15 characters or less) - * pname - the name of the parameter (25 characters or less) - * - **************************************************************************/ -void checkcell_1d(const mxArray *param, char *fname, char *pname); - - -#endif -
--- a/Problems/private/private/printf.m Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,26 +0,0 @@ -function str = printf(varargin) -%PRINTF Print formatted text to screen. -% PRINTF(FMT,VAL1,VAL2,...) formats the data in VAL1,VAL2,... according to -% the format string FMT, and prints the result to the screen. -% -% The call to PRINTF(FMT,VAL1,VAL2,...) simply invokes the call -% DISP(SPRINTF(FMT,VAL1,VAL2,...)). For a complete description of the -% format string options see function SPRINTF. -% -% STR = PRINTF(...) also returns the formatted string. - - -% Ron Rubinstein -% Computer Science Department -% Technion, Haifa 32000 Israel -% ronrubin@cs -% -% April 2008 - - -if (nargout>0) - str = sprintf(varargin{:}); - disp(str); -else - disp(sprintf(varargin{:})); -end
--- a/Problems/private/pwsmoothfield.m Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,32 +0,0 @@ -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 -
--- a/Problems/private/reggrid.m Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,136 +0,0 @@ -function [varargout] = reggrid(sz,num,mode) -%REGGRID Regular sampling grid. -% [I1,I2,...,Ip] = REGGRID([N1 N2 ... Np], NUM) returns the indices -% of a regular uniform sampling grid over a p-dimensional matrix with -% dimensions N1xN2x...xNp. NUM is the minimal number of required samples, -% and it is ensured that the actual number of samples, given by -% length(I1)xlength(I2)x...xlength(Ip), is at least as large as NUM. -% -% [I1,I2,...,Ip] = REGGRID([N1 N2 ... Np], NUM,'MODE') specifies the -% method for distributing the samples along each dimension. Valid modes -% include 'eqdist' (the default mode) and 'eqnum'. 'eqdist' indicates an -% equal distance between the samples in each dimension, while 'eqnum' -% indicates an equal number of samples in each dimension. -% -% Notes about MODE: -% -% 1. The 'eqnum' mode will generally fail when the p-th root of NUM -% (i.e. NUM^(1/p)) is larger than min([N1 N2 ... Np]). Thus 'eqdist' is -% the more useful choice for sampling an arbitrary number of samples -% from the matrix (up to the total number of matrix entries). -% -% 2. In both modes, the equality (of the distance between samples, or -% the number of samples in each dimension) is only approximate. This is -% because REGGRID attempts to maintain the appropriate equality while at -% the same time find a sampling pattern where the total number of -% samples is as close as possible to NUM. In general, the larger {Ni} -% and NUM are, the tighter the equality. -% -% Example: Sample a set of blocks uniformly from a 2D image. -% -% n = 512; blocknum = 20000; blocksize = [8 8]; -% im = rand(n,n); -% [i1,i2] = reggrid(size(im)-blocksize+1, blocknum); -% blocks = sampgrid(im, blocksize, i1, i2); -% -% See also SAMPGRID. - -% Ron Rubinstein -% Computer Science Department -% Technion, Haifa 32000 Israel -% ronrubin@cs -% -% November 2007 - -dim = length(sz); - -if (nargin<3) - mode = 'eqdist'; -end - -if (any(sz<1)) - error(['Invalid matrix size : [' num2str(sz) ']']); -end - -if (num > prod(sz)) - warning(['Invalid number of samples, returning maximum number of samples.']); -elseif (num <= 0) - if (num < 0) - warning('Invalid number of samples, assuming 0 samples.'); - end - for i = 1:length(sz) - varargout{i} = []; - end - return; -end - - -if (strcmp(mode,'eqdist')) - - % approximate distance between samples: total volume divided by number of - % samples gives the average volume per sample. then, taking the p-th root - % gives the average distance between samples - d = (prod(sz)/num)^(1/dim); - - % compute the initial guess for number of samples in each dimension. - % then, while total number of samples is too large, decrese the number of - % samples by one in the dimension where the samples are the most crowded. - % finally, do the opposite process until just passing num, so the final - % number of samples is the closest to num from above. - - n = min(max(round(sz/d),1),sz); % set n so that it saturates at 1 and sz - - active_dims = find(n>1); % dimensions where the sample num can be reduced - while(prod(n)>num && ~isempty(active_dims)) - [y,id] = min((sz(active_dims)-1)./n(active_dims)); - n(active_dims(id)) = n(active_dims(id))-1; - if (n(active_dims(id)) < 2) - active_dims = find(n>1); - end - end - - active_dims = find(n<sz); % dimensions where the sample num can be increased - while(prod(n)<num && ~isempty(active_dims)) - [y,id] = max((sz(active_dims)-1)./n(active_dims)); - n(active_dims(id)) = n(active_dims(id))+1; - if (n(active_dims(id)) >= sz(active_dims(id))) - active_dims = find(n<sz); - end - end - - for i = 1:dim - varargout{i} = round((1:n(i))/n(i)*sz(i)); - varargout{i} = varargout{i} - floor((varargout{i}(1)-1)/2); - end - -elseif (strcmp(mode,'eqnum')) - - % same idea as above - n = min(max( ones(size(sz)) * round(num^(1/dim)) ,1),sz); - - active_dims = find(n>1); - while(prod(n)>num && ~isempty(active_dims)) - [y,id] = min((sz(active_dims)-1)./n(active_dims)); - n(active_dims(id)) = n(active_dims(id))-1; - if (n(active_dims(id)) < 2) - active_dims = find(n>1); - end - end - - active_dims = find(n<sz); - while(prod(n)<num && ~isempty(active_dims)) - [y,id] = max((sz(active_dims)-1)./n(active_dims)); - n(active_dims(id)) = n(active_dims(id))+1; - if (n(active_dims(id)) >= sz(active_dims(id))) - active_dims = find(n<sz); - end - end - - for i = 1:dim - varargout{i} = round((1:n(i))/n(i)*sz(i)); - varargout{i} = varargout{i} - floor((varargout{i}(1)-1)/2); - end -else - error('Invalid sampling mode'); -end -
--- a/Problems/private/remove_dc.m Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,31 +0,0 @@ -function [y,dc] = remove_dc(x,columns) -%REMOVE_DC Remove DC channel from signals. -% [Y,DC] = REMOVE_DC(X) removes the DC channel (i.e. the mean) from the -% specified (possibly multi-dimensional) signal X. Y is the DC-free -% signal and is the same size as X. DC is a scalar containing the mean of -% the signal X. -% -% [Y,DC] = REMOVE_DC(X,'columns') where X is a 2D matrix, treats the -% columns of X as a set of 1D signals, removing the DC channel from each -% one individually. Y is the same size as X and contains the DC-free -% signals. DC is a row vector of length size(X,2) containing the means of -% the signals in X. -% -% See also ADD_DC. - - -if (nargin==2 && strcmpi(columns,'columns')), columns = 1; -else columns = 0; -end - -if (columns) - dc = mean(x); - y = addtocols(x,-dc); -else - if (ndims(x)==2) % temporary, will remove in future - warning('Treating 2D matrix X as a single signal and not each column individually'); - end - dc = mean(x(:)); - y = x-dc; -end -
--- a/Problems/private/rowlincomb.c Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,148 +0,0 @@ -/************************************************************************** - * - * File name: rowlincomb.c - * - * Ron Rubinstein - * Computer Science Department - * Technion, Haifa 32000 Israel - * ronrubin@cs - * - * Last Updated: 21.5.2009 - * - *************************************************************************/ - -#include "mex.h" - - -/* Input Arguments */ - -#define X_IN prhs[0] -#define A_IN prhs[1] -#define ROWS_IN prhs[2] -#define COLS_IN prhs[3] - - -/* Output Arguments */ - -#define Y_OUT plhs[0] - - -void mexFunction(int nlhs, mxArray *plhs[], -int nrhs, const mxArray*prhs[]) - -{ - double *A, *x, *y, *rows, *cols; - mwSize m,n,m1,n1,m2,n2,rownum,colnum; - mwIndex *row_ids,*col_ids,i,j; - int colnumspecified=0; - - - /* Check for proper number of arguments */ - - if (nrhs!=3 && nrhs!=4) { - mexErrMsgTxt("Invalid number of input arguments."); - } else if (nlhs > 1) { - mexErrMsgTxt("Too many output arguments."); - } - - - /* Check the input dimensions */ - - m = mxGetM(A_IN); - n = mxGetN(A_IN); - if (!mxIsDouble(A_IN) || mxIsComplex(A_IN) || mxGetNumberOfDimensions(A_IN)>2) { - mexErrMsgTxt("ROWLINCOMB requires that A be a double matrix."); - } - - m1 = mxGetM(ROWS_IN); - n1 = mxGetN(ROWS_IN); - if (!mxIsDouble(ROWS_IN) || mxIsComplex(ROWS_IN) || (m1!=1 && n1!=1)) { - mexErrMsgTxt("ROWLINCOMB requires that ROWS be an index vector of type double."); - } - rownum = (m1 > n1) ? m1 : n1; /* the number of rows in the linear combination */ - - m2 = mxGetM(X_IN); - n2 = mxGetN(X_IN); - if (!mxIsDouble(X_IN) || mxIsComplex(X_IN) || ((m2!=1) && (n2!=1))) { - mexErrMsgTxt("ROWLINCOMB requires that X be a double vector."); - } - - if (m2 != rownum && n2 != rownum) { - mexErrMsgTxt("The length of X does not match the number of rows in ROWS."); - } - - if (nrhs==4) { - m1 = mxGetM(COLS_IN); - n1 = mxGetN(COLS_IN); - if (!mxIsDouble(COLS_IN) || mxIsComplex(COLS_IN) || (m1!=1 && n1!=1)) { - mexErrMsgTxt("ROWLINCOMB requires that COLS be an index vector of type double."); - } - colnum = (m1 > n1) ? m1 : n1; /* the number of columns */ - colnumspecified = 1; - cols = mxGetPr(COLS_IN); - - Y_OUT = mxCreateDoubleMatrix(1, colnum, mxREAL); - } - else { - cols = 0; - Y_OUT = mxCreateDoubleMatrix(1, n, mxREAL); - } - - - /* Assign pointers to the various parameters */ - A = mxGetPr(A_IN); - rows = mxGetPr(ROWS_IN); - x = mxGetPr(X_IN); - y = mxGetPr(Y_OUT); - - - /* check row indices */ - - row_ids = (mwIndex*)mxMalloc(rownum*sizeof(mwIndex)); - - for (i=0; i<rownum; ++i) { - row_ids[i] = (mwIndex)(rows[i]+0.1)-1; - if (row_ids[i]<0 || row_ids[i]>=m) { - mexErrMsgTxt("Row index in ROWS is out of range."); - } - } - - - - if (colnumspecified) { - - /* check column indices */ - col_ids = (mwIndex*)mxMalloc(colnum*sizeof(mwIndex)); - - for (i=0; i<colnum; ++i) { - col_ids[i] = (mwIndex)(cols[i]+0.1)-1; - if (col_ids[i]<0 || col_ids[i]>=n) { - mexErrMsgTxt("Column index in COLS is out of range."); - } - } - - /* Do the actual computation */ - for (j=0; j<colnum; ++j) { - for (i=0; i<rownum; ++i) { - y[j] += A[m*col_ids[j]+row_ids[i]]*x[i]; - } - } - - mxFree(col_ids); - } - - else { - - /* Do the actual computation */ - for (j=0; j<n; ++j) { - for (i=0; i<rownum; ++i) { - y[j] += A[m*j+row_ids[i]]*x[i]; - } - } - } - - - mxFree(row_ids); - - return; -}
--- a/Problems/private/rowlincomb.m Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,26 +0,0 @@ -%ROWLINCOMB Linear combination of matrix rows. -% Y = ROWLINCOMB(X,A,ROWS) computes a linear combination of the rows of -% the matrix A. The row indices are specified in the vector ROWS, and the -% correspoinding coefficients are specified in the vector X. The vectors -% ROWS and X must be of the same length. The call Y = ROWLINCOMB(X,A,ROWS) -% is essentially equivalent to the command -% -% Y = X'*A(ROWS,:) . -% -% However, it is implemented much more efficiently. -% -% Y = ROWLINCOMB(X,A,ROWS,COLS) only works on the columns of A specified -% in COLS, returning a vector of length equal to COLS. This call is -% essentially equivalent to the command -% -% Y = X'*A(ROWS,COLS) . -% -% See also COLLINCOMB. - - -% Ron Rubinstein -% Computer Science Department -% Technion, Haifa 32000 Israel -% ronrubin@cs -% -% April 2009
--- a/Problems/private/sampgrid.m Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,72 +0,0 @@ -function y = sampgrid(x,blocksize,varargin) -%SAMPGRID Sample a multi-dimensional matrix on a regular grid. -% Y = SAMPGRID(X,BLOCKSIZE,I1,I2,...,Ip) extracts block samples of size -% BLOCKSIZE from the p-dimensional matrix X, arranging the samples as the -% column vectors of the matrix Y. The locations of the (1,1,..,1)-th -% elements of each block are given in the index vectors I1,I2,..Ip. The -% total number of samples taken is length(I1)xlength(I2)x...xlength(Ip). -% BLOCKSIZE should either be a p-element vector of the form [N1,N2,...Np], -% or a scalar N which is shorthand for the square block size [N N ... N]. -% -% Example: Sample a set of blocks uniformly from a 2D image. -% -% n = 512; blocknum = 20000; blocksize = [8 8]; -% im = rand(n,n); -% [i1,i2] = reggrid(size(im)-blocksize+1, blocknum); -% blocks = sampgrid(im, blocksize, i1, i2); -% -% See also REGGRID. - -% Ron Rubinstein -% Computer Science Department -% Technion, Haifa 32000 Israel -% ronrubin@cs -% -% November 2007 - - -p = ndims(x); -if (p==2 && any(size(x)==1) && length(blocksize)==1) - p = 1; -end - -if (numel(blocksize)==1) - blocksize = ones(1,p)*blocksize; -end - -n = zeros(1,p); -for i = 1:p - n(i) = length(varargin{i}); -end - -nsamps = prod(n); - -% create y of the same class as x -y = zeros(prod(blocksize),nsamps,class(x)); - -% ids() contains the index of the current block in I1..Ip -ids = ones(p,1); - -% block_ids contains the indices of the current block in X -block_ids = cell(p,1); -for j = 1:p - block_ids{j} = varargin{j}(1) : varargin{j}(1)+blocksize(j)-1; -end - -for k = 1:nsamps - block = x(block_ids{:}); - y(:,k) = block(:); - - % increment ids() and block_ids{} - if (k<nsamps) - j = 1; - while (ids(j) == n(j)) - ids(j) = 1; - block_ids{j} = varargin{j}(1) : varargin{j}(1)+blocksize(j)-1; - j = j+1; - end - ids(j) = ids(j)+1; - block_ids{j} = varargin{j}(ids(j)) : varargin{j}(ids(j))+blocksize(j)-1; - end -end -
--- a/Problems/private/scalarToRGB.m Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,18 +0,0 @@ -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
--- a/Problems/private/secs2hms.m Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,28 +0,0 @@ -function [h,m,s] = secs2hms(t) -%SECS2HMS Convert seconds to hours, minutes and seconds. -% [H,M,S] = SECS2HMS(T) converts the specified number of seconds T to -% hours, minutes and seconds. H and M are whole numbers, and S is real. -% -% Example: Estimate the remaining time of a loop -% -% n = 10; tic; -% for i = 1:n -% pause(1); -% [h,m,s] = secs2hms( (n-i)*toc/i ); -% printf('estimated remaining time: %02d:%02d:%05.2f',h,m,s); -% end - - -% Ron Rubinstein -% Computer Science Department -% Technion, Haifa 32000 Israel -% ronrubin@cs -% -% April 2008 - - -s = t; -h = fix(s/3600); -s = rem(s,3600); -m = fix(s/60); -s = rem(s,60);
--- a/Problems/private/spdiag.m Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,38 +0,0 @@ -function Y = spdiag(V,K) -%SPDIAG Sparse diagonal matrices. -% SPDIAG(V,K) when V is a vector with N components is a sparse square -% matrix of order N+ABS(K) with the elements of V on the K-th diagonal. -% K = 0 is the main diagonal, K > 0 is above the main diagonal and K < 0 -% is below the main diagonal. -% -% SPDIAG(V) is the same as SPDIAG(V,0) and puts V on the main diagonal. -% -% See also DIAG, SPDIAGS. - - -% Ron Rubinstein -% Computer Science Department -% Technion, Haifa 32000 Israel -% ronrubin@cs -% -% June 2008 - - -if (nargin<2) - K = 0; -end - -n = length(V) + abs(K); - -if (K>0) - i = 1:length(V); - j = K+1:n; -elseif (K<0) - i = -K+1:n; - j = 1:length(V); -else - i = 1:n; - j = 1:n; -end - -Y = sparse(i,j,V(:),n,n); \ No newline at end of file
--- a/Problems/private/sprow.c Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,94 +0,0 @@ -/************************************************************************** - * - * File name: sprow.c - * - * Ron Rubinstein - * Computer Science Department - * Technion, Haifa 32000 Israel - * ronrubin@cs - * - * Last Updated: 24.8.2009 - * - *************************************************************************/ - - -#include "mex.h" -#include "mexutils.h" - - -/* Input Arguments */ - -#define A_IN prhs[0] -#define J_IN prhs[1] - - -/* Output Arguments */ - -#define X_OUT plhs[0] -#define ID_OUT plhs[1] - - -void mexFunction(int nlhs, mxArray *plhs[], - int nrhs, const mxArray*prhs[]) - -{ - double *pr, *x, *id, rowid; - mwIndex *ir, *jc; - mwSize m, n; - mwIndex i, j, k, l, rowlen; - - if (nrhs != 2) { - mexErrMsgTxt("GETSPROW requires two input arguments."); - } else if (nlhs > 2) { - mexErrMsgTxt("Too many output arguments."); - } - - checkmatrix(A_IN, "GETSPROW", "A"); - checksparse(A_IN, "GETSPROW", "A"); - checkscalar(J_IN, "GETSPROW", "J"); - - m = mxGetM(A_IN); - n = mxGetN(A_IN); - - rowid = mxGetScalar(J_IN); - if (rowid < 0) { - mexErrMsgTxt("Invalid row index."); - } - j = (mwIndex)(rowid + 1e-2); - if (j<1 || j>m) { - mexErrMsgTxt("Row index is out of range."); - } - j--; - - pr = mxGetPr(A_IN); - ir = mxGetIr(A_IN); - jc = mxGetJc(A_IN); - - /* Determine length of row */ - rowlen = 0; - for (i=0; i<jc[n]; ++i) { - rowlen += (ir[i]==j); - } - - /* Allocate output parameters */ - X_OUT = mxCreateDoubleMatrix(1, rowlen, mxREAL); - ID_OUT = mxCreateDoubleMatrix(1, rowlen, mxREAL); - - x = mxGetPr(X_OUT); - id = mxGetPr(ID_OUT); - - /* Compute j-th row */ - k=0; - for (l=1; l<=n; ++l) { - i = jc[l-1]; - while (i<jc[l] && ir[i]<j) { - i++; - } - if (i<jc[l] && ir[i]==j) { - x[k] = pr[i]; - id[k] = l; - k++; - } - } - -}
--- a/Problems/private/sprow.m Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,19 +0,0 @@ -%SPROW Extract row of sparse matrix. -% X = SPROW(A,J) where A is a sparse matrix, returns the nonzero values in -% the row A(J,:). -% -% [X,ID] = SPROW(A,J) also returns the column indices of the nonzeros. -% -% Note that the call [X,ID] = SPROW(A,J) is equivalent (but more efficient -% than) the Matlab code -% -% IDS = find(A(J,:)); -% X = A(J,IDS); - - -% Ron Rubinstein -% Computer Science Department -% Technion, Haifa 32000 Israel -% ronrubin@cs -% -% August 2009
--- a/Problems/private/thumbFromOp.m Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,42 +0,0 @@ -function P = thumbFromOp(op,m,n,sm,sn,grayscale) -% Output matrix P, of size m x n, sampled at -% sm x sn upper top - -% Copyright 2008, Ewout van den Berg and Michael P. Friedlander -% http://www.cs.ubc.ca/labs/scl/sparco -% $Id: thumbFromOp.m 1040 2008-06-26 20:29:02Z ewout78 $ - -if nargin < 6, grayscale = 0; end - -info = op([],0); - -sm = min(info{1},sm); -sn = min(info{2},sn); -M = zeros(sm,sn); - -for i=1:sn - v = zeros(info{2},1); v(i) = 1; - w = real(op(v,1)); - - M(:,i) = w(1:sm); -end - -mn = min(min(M)); -mx = max(max(M)); -M = (M - mn) / (mx-mn); - -idxm = floor(linspace(1,sm+1,m+1)); idxm = idxm(1:end-1); -idxn = floor(linspace(1,sn+1,n+1)); idxn = idxn(1:end-1); - -if grayscale - P = 1-M(idxm,idxn); -else - clrmap = hsv; - M = 1 + round(M * (length(clrmap)-1)); - P = zeros(m,n,3); - for j1=1:m - for j2=1:n - P(j1,j2,:) = clrmap(M(idxm(j1),idxn(j2)),:); - end - end -end
--- a/Problems/private/thumbPlot.m Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,38 +0,0 @@ -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
--- a/Problems/private/thumbwrite.m Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,18 +0,0 @@ -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);
--- a/Problems/private/timerclear.m Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,37 +0,0 @@ -function timerclear() -%TIMERCLEAR Clear all timers. -% TIMERCLEAR clears all currenly registered timers, invalidating all -% timer ids. -% -% Note: since registered timers do not consume CPU power except for when -% the TIMER<*> functions are called, this function is only useful in -% situations where a large number of timers have been initialized, and -% there is a need to reclaim memory. -% -% See also TIMERINIT, TIMERETA. - - -% Ron Rubinstein -% Computer Science Department -% Technion, Haifa 32000 Israel -% ronrubin@cs -% -% June 2008 - - -global utiltbx_timer_start_times % start times -global utiltbx_time_lastdisp % last display times -global utiltbx_timer_iternums % iteration numbers -global utiltbx_timer_lastiter % last queried iteration numbers -global utiltbx_timer_name % timer names -global utiltbx_timer_callfun % timer calling functions - - -% clear all timers % - -utiltbx_timer_start_times = []; -utiltbx_time_lastdisp = []; -utiltbx_timer_iternums = []; -utiltbx_timer_lastiter = []; -utiltbx_timer_name = []; -utiltbx_timer_callfun = [];
--- a/Problems/private/timereta.m Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,98 +0,0 @@ -function varargout = timereta(tid,iter,delay) -%TIMERETA Estimated remaining time. -% S = TIMERETA(TID,ITER) returns the estimated remaining time (in -% seconds) for the process associated with timer TID, assuming the -% process has completed ITER iterations. Note: the function will exit -% with an error if the timer TID does not have an associated number of -% iterations (see function TIMERINIT). -% -% [H,M,S] = TIMERETA(TID,ITER) returns the estimated remaining time in -% hours, minutes and seconds. -% -% TIMERETA(TID,ITER), with no output arguments, prints the estimated -% remaining time to the screen. The time is displayed in the format -% -% TIMERNAME: iteration ITER / ITERNUM, estimated remaining time: HH:MM:SS.SS -% -% If the timer has no assigned name, the display format changes to -% -% Iteration ITER / ITERNUM, estimated remaining time: HH:MM:SS.SS -% -% TIMERETA(TID,ITER,DELAY) only displays the remaining time if the -% time elapsed since the previous printout is at least DELAY seconds. -% -% See also TIMERINIT, TIMERCLEAR. - - -% Ron Rubinstein -% Computer Science Department -% Technion, Haifa 32000 Israel -% ronrubin@cs -% -% June 2008 - - -global utiltbx_timer_start_times -global utiltbx_timer_iternums -global utiltbx_timer_lastiter -global utiltbx_time_lastdisp -global utiltbx_timer_name - - -if (tid<1 || tid>length(utiltbx_timer_iternums)) - error('Unknown timer id'); -end - -if (utiltbx_timer_iternums(tid) < 0) - error('Specified timer does not have an associated number of iterations'); -end - -% update last reported iteration number -utiltbx_timer_lastiter(tid) = iter; - -% compute elapsed time -starttime = utiltbx_timer_start_times(tid,:); -currtime = clock; -timediff = etime(currtime, starttime); - -% total iteration number -iternum = utiltbx_timer_iternums(tid); - -% compute eta -timeremain = (iternum-iter)*timediff/iter; - -% return eta in seconds -if (nargout==1) - varargout{1} = timeremain; - -% return eta in hms -elseif (nargout==3) - [varargout{1}, varargout{2}, varargout{3}] = secs2hms(timeremain); - - -% print eta -elseif (nargout==0) - - % check last display time - lastdisptime = utiltbx_time_lastdisp(tid,:); - if (nargin>2 && etime(currtime,lastdisptime) < delay) - return; - end - - % update last display time - utiltbx_time_lastdisp(tid,:) = currtime; - - % display timer - [hrs,mins,secs] = secs2hms(timeremain); - if (isempty(utiltbx_timer_name{tid})) - printf('Iteration %d / %d, estimated remaining time: %02d:%02d:%05.2f', iter, iternum, hrs, mins, secs); - else - timername = utiltbx_timer_name{tid}; - printf('%s: iteration %d / %d, estimated remaining time: %02d:%02d:%05.2f', timername, iter, iternum, hrs, mins, secs); - end - -% invalid number of outputs -else - error('Invalid number of output arguments'); -end -
--- a/Problems/private/timerinit.m Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,110 +0,0 @@ -function tid = timerinit(par1,par2) -%TIMERINIT Initialize a new timer. -% TID = TIMERINIT() initializes a new timer for counting elapsed time, -% and returns its id. -% -% TID = TIMERINIT('TIMERNAME') sets the timer name to the specified -% string for display purposes. -% -% TID = TIMERINIT(ITERNUM) initializes a new ETA timer for a process with -% ITERNUM iterations. An ETA timer can be used for both counting elapsed -% time and estimating remaining time. -% -% TID = TIMERINIT('TIMERNAME',ITERNUM) sets the ETA timer name to the -% specified string for display purposes. -% -% Example: -% -% tid = timerinit(100); -% for i = 1:100 -% pause(0.07); -% timereta(tid,i,1); -% end -% timereta(tid,i); -% -% See also TIMERETA, TIMERCLEAR. - - -% Ron Rubinstein -% Computer Science Department -% Technion, Haifa 32000 Israel -% ronrubin@cs -% -% June 2008 - - -global utiltbx_timer_start_times % start times -global utiltbx_time_lastdisp % last display times -global utiltbx_timer_iternums % iteration numbers -global utiltbx_timer_lastiter % last queried iteration numbers -global utiltbx_timer_name % timer names -global utiltbx_timer_callfun % timer calling functions - - -% parse function arguments % - -if (nargin==0) - - iternum = -1; - timername = ''; - -elseif (nargin==1) - - if (ischar(par1)) - iternum = -1; - timername = par1; - - elseif (isnumeric(par1) && numel(par1)==1 && par1>0) - iternum = par1; - timername = ''; - - else - error('Invalid number of iterations'); - end - -elseif (nargin==2) - - if (ischar(par1) && isnumeric(par2)) - if (numel(par2)==1 && par2>0) - timername = par1; - iternum = par2; - else - error('Invalid number of iterations'); - end - else - error('Invalid function syntax'); - end - -else - error('Too many function parameters'); -end - - -% register the timer % - -if (isempty(utiltbx_timer_start_times)) - utiltbx_timer_start_times = clock; - utiltbx_time_lastdisp = utiltbx_timer_start_times; - utiltbx_timer_iternums = double(iternum); - utiltbx_timer_lastiter = 0; - utiltbx_timer_name = { timername }; - utiltbx_timer_callfun = {}; - tid = 1; -else - utiltbx_timer_start_times(end+1,:) = clock; - utiltbx_time_lastdisp(end+1,:) = utiltbx_timer_start_times(end,:); - utiltbx_timer_iternums(end+1) = double(iternum); - utiltbx_timer_lastiter(end+1) = 0; - utiltbx_timer_name{end+1} = timername; - tid = size(utiltbx_timer_start_times,1); -end - - -% detect timer calling function % - -st = dbstack; -if (length(dbstack) >= 2) - utiltbx_timer_callfun{end+1} = st(2).name; -else - utiltbx_timer_callfun{end+1} = ''; -end
--- a/Problems/private/updateFigure.m Mon Mar 14 17:06:56 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,93 +0,0 @@ -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