# HG changeset patch # User idamnjanovic # Date 1300123538 0 # Node ID 23f9dd7b9d785fdba6ffdde011c771d9e23620cf # Parent 66fe897f1bf4abc9b25a895b92e70f40fc9376da diff -r 66fe897f1bf4 -r 23f9dd7b9d78 Problems/private/add_dc.m --- 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 - - - diff -r 66fe897f1bf4 -r 23f9dd7b9d78 Problems/private/addtocols.c --- 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 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] 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=m) { - mexErrMsgTxt("Row index in ROWS is out of range."); - } - } - - /* Do the actual computation */ - for (i=0; i=n) { - mexErrMsgTxt("Column index in COLS is out of range."); - } - for (j=0; j=n) { - mexErrMsgTxt("Column index in COLS is out of range."); - } - for (j=0; j 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; diff -r 66fe897f1bf4 -r 23f9dd7b9d78 Problems/private/countcover.m --- 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 diff -r 66fe897f1bf4 -r 23f9dd7b9d78 Problems/private/diag_ids.m --- 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 diff -r 66fe897f1bf4 -r 23f9dd7b9d78 Problems/private/dictdist.m --- 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; diff -r 66fe897f1bf4 -r 23f9dd7b9d78 Problems/private/genPDF.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 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 - - - - - - diff -r 66fe897f1bf4 -r 23f9dd7b9d78 Problems/private/genSampling.m --- 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)) - - -/* 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]=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 diff -r 66fe897f1bf4 -r 23f9dd7b9d78 Problems/private/imnormalize.m --- 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); diff -r 66fe897f1bf4 -r 23f9dd7b9d78 Problems/private/iswhole.m --- 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; diff -r 66fe897f1bf4 -r 23f9dd7b9d78 Problems/private/make.m --- 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 diff -r 66fe897f1bf4 -r 23f9dd7b9d78 Problems/private/mexutils.c --- 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 - - - -/* 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); - } -} - diff -r 66fe897f1bf4 -r 23f9dd7b9d78 Problems/private/mexutils.h --- 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 - diff -r 66fe897f1bf4 -r 23f9dd7b9d78 Problems/private/myblas.c --- 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 - - -/* find maximum of absolute values */ - -mwIndex maxabs(double c[], mwSize m) -{ - mwIndex maxid=0, k; - double absval, maxval = SQR(*c); /* use square which is quicker than absolute value */ - - for (k=1; k maxval) { - maxval = absval; - maxid = k; - } - } - return maxid; -} - - -/* compute y := alpha*x + y */ - -void vec_sum(double alpha, double x[], double y[], mwSize n) -{ - mwIndex i; - - for (i=0; i maxval) { - maxval = val; - maxid = k; - } - } - return maxid; -} - - -/* solve L*x = b */ - -void backsubst_L(double L[], double b[], double x[], mwSize n, mwSize k) -{ - mwIndex i, j; - double rhs; - - for (i=0; i=1; --i) { - rhs = b[i-1]; - for (j=i; j=1; --i) { - rhs = b[i-1]; - for (j=i; j - - - -/************************************************************************** - * Squared value. - **************************************************************************/ -#define SQR(X) ((X)*(X)) - - - -/************************************************************************** - * Matrix-vector multiplication. - * - * Computes an operation of the form: - * - * y := alpha*A*x - * - * Parameters: - * A - matrix of size n X m - * x - vector of length m - * y - output vector of length n - * alpha - real constant - * n, m - dimensions of A - * - * Note: This function re-writes the contents of y. - * - **************************************************************************/ -void mat_vec(double alpha, double A[], double x[], double y[], mwSize n, mwSize m); - - - -/************************************************************************** - * Matrix-transpose-vector multiplication. - * - * Computes an operation of the form: - * - * y := alpha*A'*x - * - * Parameters: - * A - matrix of size n X m - * x - vector of length n - * y - output vector of length m - * alpha - real constant - * n, m - dimensions of A - * - * Note: This function re-writes the contents of y. - * - **************************************************************************/ -void matT_vec(double alpha, double A[], double x[], double y[], mwSize n, mwSize m); - - - -/************************************************************************** - * Sparse-matrix-vector multiplication. - * - * Computes an operation of the form: - * - * y := alpha*A*x - * - * where A is a sparse matrix. - * - * Parameters: - * pr,ir,jc - sparse representation of the matrix A, of size n x m - * x - vector of length m - * y - output vector of length n - * alpha - real constant - * n, m - dimensions of A - * - * Note: This function re-writes the contents of y. - * - **************************************************************************/ -void mat_sp_vec(double alpha, double pr[], mwIndex ir[], mwIndex jc[], double x[], double y[], mwSize n, mwSize m); - - - -/************************************************************************** - * Sparse-matrix-transpose-vector multiplication. - * - * Computes an operation of the form: - * - * y := alpha*A'*x - * - * where A is a sparse matrix. - * - * Parameters: - * pr,ir,jc - sparse representation of the matrix A, of size n x m - * x - vector of length m - * y - output vector of length n - * alpha - real constant - * n, m - dimensions of A - * - * Note: This function re-writes the contents of y. - * - **************************************************************************/ -void matT_sp_vec(double alpha, double pr[], mwIndex ir[], mwIndex jc[], double x[], double y[], mwSize n, mwSize m); - - - -/************************************************************************** - * Matrix-sparse-vector multiplication. - * - * Computes an operation of the form: - * - * y := alpha*A*x - * - * where A is a matrix and x is a sparse vector. - * - * Parameters: - * A - matrix of size n X m - * pr,ir,jc - sparse representation of the vector x, of length m - * y - output vector of length n - * alpha - real constant - * n, m - dimensions of A - * - * Note: This function re-writes the contents of y. - * - **************************************************************************/ -void mat_vec_sp(double alpha, double A[], double pr[], mwIndex ir[], mwIndex jc[], double y[], mwSize n, mwSize m); - - - -/************************************************************************** - * Matrix-transpose-sparse-vector multiplication. - * - * Computes an operation of the form: - * - * y := alpha*A'*x - * - * where A is a matrix and x is a sparse vector. - * - * Parameters: - * A - matrix of size n X m - * pr,ir,jc - sparse representation of the vector x, of length n - * y - output vector of length m - * alpha - real constant - * n, m - dimensions of A - * - * Note: This function re-writes the contents of y. - * - **************************************************************************/ -void matT_vec_sp(double alpha, double A[], double pr[], mwIndex ir[], mwIndex jc[], double y[], mwSize n, mwSize m); - - - -/************************************************************************** - * Sparse-matrix-sparse-vector multiplication. - * - * Computes an operation of the form: - * - * y := alpha*A*x - * - * where A is a sparse matrix and x is a sparse vector. - * - * Parameters: - * pr,ir,jc - sparse representation of the matrix A, of size n x m - * prx,irx,jcx - sparse representation of the vector x (of length m) - * y - output vector of length n - * alpha - real constant - * n, m - dimensions of A - * - * Note: This function re-writes the contents of y. - * - **************************************************************************/ -void mat_sp_vec_sp(double alpha, double pr[], mwIndex ir[], mwIndex jc[], double prx[], mwIndex irx[], mwIndex jcx[], double y[], mwSize n, mwSize m); - - - -/************************************************************************** - * Sparse-matrix-transpose-sparse-vector multiplication. - * - * Computes an operation of the form: - * - * y := alpha*A'*x - * - * where A is a sparse matrix and x is a sparse vector. - * - * Importnant note: this function is provided for completeness, but is NOT efficient. - * If possible, convert x to non-sparse representation and use matT_vec_sp instead. - * - * Parameters: - * pr,ir,jc - sparse representation of the matrix A, of size n x m - * prx,irx,jcx - sparse representation of the vector x (of length n) - * y - output vector of length n - * alpha - real constant - * n, m - dimensions of A - * - * Note: This function re-writes the contents of y. - * - **************************************************************************/ -void matT_sp_vec_sp(double alpha, double pr[], mwIndex ir[], mwIndex jc[], double prx[], mwIndex irx[], mwIndex jcx[], double y[], mwSize n, mwSize m); - - - -/************************************************************************** - * Matrix-matrix multiplication. - * - * Computes an operation of the form: - * - * X := alpha*A*B - * - * Parameters: - * A - matrix of size n X m - * B - matrix of size m X k - * X - output matrix of size n X k - * alpha - real constant - * n, m, k - dimensions of A, B - * - * Note: This function re-writes the contents of X. - * - **************************************************************************/ -void mat_mat(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k); - - - -/************************************************************************** - * Matrix-transpose-matrix multiplication. - * - * Computes an operation of the form: - * - * X := alpha*A*B - * - * Parameters: - * A - matrix of size n X m - * B - matrix of size m X k - * X - output matrix of size n X k - * alpha - real constant - * n, m, k - dimensions of A, B - * - * Note: This function re-writes the contents of X. - * - **************************************************************************/ -void matT_mat(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k); - - - -/************************************************************************** - * Tensor-matrix multiplication. - * - * This function accepts a 3-D tensor A of size n X m X k - * and a 2-D matrix B of size l X k. - * The function computes the 3-D tensor X of size n X m X l, where - * - * X(i,j,:) = B*A(i,j,:) - * - * for all i,j. - * - * Parameters: - * A - tensor of size n X m X k - * B - matrix of size l X k - * X - output tensor of size n X m X l - * alpha - real constant - * n, m, k, l - dimensions of A, B - * - * Note: This function re-writes the contents of X. - * - **************************************************************************/ -void tens_mat(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k, mwSize l); - - - -/************************************************************************** - * Tensor-matrix-transpose multiplication. - * - * This function accepts a 3-D tensor A of size n X m X k - * and a 2-D matrix B of size k X l. - * The function computes the 3-D tensor X of size n X m X l, where - * - * X(i,j,:) = B'*A(i,j,:) - * - * for all i,j. - * - * Parameters: - * A - tensor of size n X m X k - * B - matrix of size k X l - * X - output tensor of size n X m X l - * alpha - real constant - * n, m, k, l - dimensions of A, B - * - * Note: This function re-writes the contents of X. - * - **************************************************************************/ -void tens_matT(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k, mwSize l); - - - -/************************************************************************** - * Vector-vector sum. - * - * Computes an operation of the form: - * - * y := alpha*x + y - * - * Parameters: - * x - vector of length n - * y - output vector of length n - * alpha - real constant - * n - length of x,y - * - * Note: This function re-writes the contents of y. - * - **************************************************************************/ -void vec_sum(double alpha, double x[], double y[], mwSize n); - - - -/************************************************************************** - * Triangular back substitution. - * - * Solve the set of linear equations - * - * T*x = b - * - * where T is lower or upper triangular. - * - * Parameters: - * ul - 'U' for upper triangular, 'L' for lower triangular - * A - matrix of size n x m containing T - * b - vector of length k - * x - output vector of length k - * n - size of first dimension of A - * k - the size of the equation set, k<=n,m - * - * Note: - * The matrix A can be of any size n X m, as long as n,m >= k. - * Only the lower/upper triangle of the submatrix A(1:k,1:k) defines the - * matrix T (depending on the parameter ul). - * - **************************************************************************/ -void backsubst(char ul, double A[], double b[], double x[], mwSize n, mwSize k); - - - -/************************************************************************** - * Solve a set of equations using a Cholesky decomposition. - * - * Solve the set of linear equations - * - * M*x = b - * - * where M is positive definite with a known Cholesky decomposition: - * either M=L*L' (L lower triangular) or M=U'*U (U upper triangular). - * - * Parameters: - * ul - 'U' for upper triangular, 'L' for lower triangular decomposition - * A - matrix of size n x m with the Cholesky decomposition of M - * b - vector of length k - * x - output vector of length k - * n - size of first dimension of A - * k - the size of the equation set, k<=n,m - * - * Note: - * The matrix A can be of any size n X m, as long as n,m >= k. - * Only the lower/upper triangle of the submatrix A(1:k,1:k) is used as - * the Cholesky decomposition of M (depending on the parameter ul). - * - **************************************************************************/ -void cholsolve(char ul, double A[], double b[], double x[], mwSize n, mwSize k); - - - -/************************************************************************** - * Maximum absolute value. - * - * Returns the index of the coefficient with maximal absolute value in a vector. - * - * Parameters: - * x - vector of length n - * n - length of x - * - **************************************************************************/ -mwIndex maxabs(double x[], mwSize n); - - - -/************************************************************************** - * Maximum vector element. - * - * Returns the index of the maximal coefficient in a vector. - * - * Parameters: - * x - vector of length n - * n - length of x - * - **************************************************************************/ -mwIndex maxpos(double x[], mwSize n); - - - -/************************************************************************** - * Vector-vector dot product. - * - * Computes an operation of the form: - * - * c = a'*b - * - * Parameters: - * a, b - vectors of length n - * n - length of a,b - * - * Returns: The dot product c. - * - **************************************************************************/ -double dotprod(double a[], double b[], mwSize n); - - - -/************************************************************************** - * Indexed vector assignment. - * - * Perform a permutation assignment of the form - * - * y = x(ind) - * - * where ind is an array of indices to x. - * - * Parameters: - * y - output vector of length k - * x - input vector of arbitrary length - * ind - array of indices into x (indices begin at 0) - * k - length of the array ind - * - **************************************************************************/ -void vec_assign(double y[], double x[], mwIndex ind[], mwSize k); - - - -/************************************************************************** - * Matrix transpose. - * - * Computes Y := X' - * - * Parameters: - * X - input matrix of size n X m - * Y - output matrix of size m X n - * n, m - dimensions of X - * - **************************************************************************/ -void transpose(double X[], double Y[], mwSize n, mwSize m); - - - -/************************************************************************** - * Print a matrix. - * - * Parameters: - * A - matrix of size n X m - * n, m - dimensions of A - * matname - name of matrix to display - * - **************************************************************************/ -void printmat(double A[], int n, int m, char* matname); - - - -/************************************************************************** - * Print a sparse matrix. - * - * Parameters: - * A - sparse matrix of type double - * matname - name of matrix to display - * - **************************************************************************/ -void printspmat(mxArray *A, char* matname); - - -#endif - diff -r 66fe897f1bf4 -r 23f9dd7b9d78 Problems/private/normcols.m --- 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))); diff -r 66fe897f1bf4 -r 23f9dd7b9d78 Problems/private/omp2mex.c --- 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; -} diff -r 66fe897f1bf4 -r 23f9dd7b9d78 Problems/private/omp2mex.m --- 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 diff -r 66fe897f1bf4 -r 23f9dd7b9d78 Problems/private/ompcore.c --- 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 -#include - - - -/****************************************************************************** - * * - * 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; signumeps2 && 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; ieps2 && 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 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= 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; j0 && (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; -} diff -r 66fe897f1bf4 -r 23f9dd7b9d78 Problems/private/ompcore.h --- 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 diff -r 66fe897f1bf4 -r 23f9dd7b9d78 Problems/private/ompmex.c --- 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; -} - diff -r 66fe897f1bf4 -r 23f9dd7b9d78 Problems/private/ompmex.m --- 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 diff -r 66fe897f1bf4 -r 23f9dd7b9d78 Problems/private/ompprof.c --- 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); -} - diff -r 66fe897f1bf4 -r 23f9dd7b9d78 Problems/private/ompprof.h --- 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 - - - -/************************************************************************** - * - * Constants and data types. - * - **************************************************************************/ - - -/* constants denoting the various parts of the algorithm */ - -enum { DtX_TIME, XtX_TIME, DtR_TIME, MAXABS_TIME, DtD_TIME, LCHOL_TIME, COMPCOEF_TIME, - UPDATE_DtR_TIME, UPDATE_RESNORM_TIME, COMPRES_TIME, INDEXSORT_TIME }; - - - -/* profiling data container with counters for each part of the algorithm */ - -typedef struct profdata -{ - clock_t prevtime; /* the time when last initialization/call to addproftime() was performed */ - - clock_t DtX_time; - clock_t XtX_time; - clock_t DtR_time; - clock_t maxabs_time; - clock_t DtD_time; - clock_t Lchol_time; - clock_t compcoef_time; - clock_t update_DtR_time; - clock_t update_resnorm_time; - clock_t compres_time; - clock_t indexsort_time; - - /* flags indicating whether profiling data was gathered */ - int DtX_time_counted; - int XtX_time_counted; - int DtR_time_counted; - int DtD_time_counted; - int update_DtR_time_counted; - int resnorm_time_counted; - int compres_time_counted; - int indexsort_time_counted; - -} profdata; - - - -/************************************************************************** - * - * Initialize a profdata structure, zero all counters, and start its timer. - * - **************************************************************************/ -void initprofdata(profdata *pd); - - -/************************************************************************** - * - * Add elapsed time from last call to addproftime(), or from initialization - * of profdata, to the counter specified by comptype. comptype must be one - * of the constants in the enumeration above. - * - **************************************************************************/ -void addproftime(profdata *pd, int comptype); - - -/************************************************************************** - * - * Print the current contents of the counters in profdata. - * - * Parameters: - * pd - the profdata to print - * erroromp - indicates whether error-based (nonzero) or sparsity-based (zero) - * omp was performed. - * batchomp - indicates whether batch-omp (nonzero) or omp-cholesky (zero) - * omp was performed. - * signum - number of signals processed by omp - * - **************************************************************************/ -void printprofinfo(profdata *pd, int erroromp, int batchomp, int signum); - - -#endif - diff -r 66fe897f1bf4 -r 23f9dd7b9d78 Problems/private/omputils.c --- 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 - - -const char FULL_GAMMA_STR[] = "full"; -const char SPARSE_GAMMA_STR[] = "sparse"; - - -/* convert seconds to hours, minutes and seconds */ - -void secs2hms(double sectot, int *hrs, int *mins, double *secs) -{ - *hrs = (int)(floor(sectot/3600)+1e-2); - sectot = sectot - 3600*(*hrs); - *mins = (int)(floor(sectot/60)+1e-2); - *secs = sectot - 60*(*mins); -} - - -/* quicksort, public-domain C implementation by Darel Rex Finley. */ -/* modification: sorts the array data[] as well, according to the values in the array vals[] */ - -#define MAX_LEVELS 300 - -void quicksort(mwIndex vals[], double data[], mwIndex n) { - - long piv, beg[MAX_LEVELS], end[MAX_LEVELS], i=0, L, R, swap ; - double datapiv; - - beg[0]=0; - end[0]=n; - - while (i>=0) { - - L=beg[i]; - R=end[i]-1; - - if (L=piv && L end[i-1]-beg[i-1]) { - swap=beg[i]; beg[i]=beg[i-1]; beg[i-1]=swap; - swap=end[i]; end[i]=end[i-1]; end[i-1]=swap; - } - } - else { - i--; - } - } -} diff -r 66fe897f1bf4 -r 23f9dd7b9d78 Problems/private/omputils.h --- 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 - diff -r 66fe897f1bf4 -r 23f9dd7b9d78 Problems/private/printf.m --- 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 diff -r 66fe897f1bf4 -r 23f9dd7b9d78 Problems/private/private/make.m --- 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{:}); - diff -r 66fe897f1bf4 -r 23f9dd7b9d78 Problems/private/private/mexutils.c --- 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 - - - -/* 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); - } -} - diff -r 66fe897f1bf4 -r 23f9dd7b9d78 Problems/private/private/mexutils.h --- 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 - diff -r 66fe897f1bf4 -r 23f9dd7b9d78 Problems/private/private/printf.m --- 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 diff -r 66fe897f1bf4 -r 23f9dd7b9d78 Problems/private/pwsmoothfield.m --- 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 - diff -r 66fe897f1bf4 -r 23f9dd7b9d78 Problems/private/reggrid.m --- 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(active_dims(id))) - active_dims = find(n1); - 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(active_dims(id))) - active_dims = find(n 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=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=n) { - mexErrMsgTxt("Column index in COLS is out of range."); - } - } - - /* Do the actual computation */ - for (j=0; j 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 diff -r 66fe897f1bf4 -r 23f9dd7b9d78 Problems/private/sprow.c --- 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 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 = []; diff -r 66fe897f1bf4 -r 23f9dd7b9d78 Problems/private/timereta.m --- 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 - diff -r 66fe897f1bf4 -r 23f9dd7b9d78 Problems/private/timerinit.m --- 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 diff -r 66fe897f1bf4 -r 23f9dd7b9d78 Problems/private/updateFigure.m --- 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