Mercurial > hg > smallbox
changeset 60:ad36f80e2ccf
(none)
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/SolveFISTA.m Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,211 @@ +% Copyright ©2010. The Regents of the University of California (Regents). +% All Rights Reserved. Contact The Office of Technology Licensing, +% UC Berkeley, 2150 Shattuck Avenue, Suite 510, Berkeley, CA 94720-1620, +% (510) 643-7201, for commercial licensing opportunities. + +% Authors: Arvind Ganesh, Allen Y. Yang, Zihan Zhou. +% Contact: Allen Y. Yang, Department of EECS, University of California, +% Berkeley. <yang@eecs.berkeley.edu> + +% IN NO EVENT SHALL REGENTS BE LIABLE TO ANY PARTY FOR DIRECT, INDIRECT, +% SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES, INCLUDING LOST PROFITS, +% ARISING OUT OF THE USE OF THIS SOFTWARE AND ITS DOCUMENTATION, EVEN IF +% REGENTS HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +% REGENTS SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT LIMITED +% TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A +% PARTICULAR PURPOSE. THE SOFTWARE AND ACCOMPANYING DOCUMENTATION, IF ANY, +% PROVIDED HEREUNDER IS PROVIDED "AS IS". REGENTS HAS NO OBLIGATION TO +% PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS. + +%% This function is modified from Matlab code proximal_gradient_bp + +function [x_hat,nIter] = SolveFISTA(A,b, varargin) + +% b - m x 1 vector of observations/data (required input) +% A - m x n measurement matrix (required input) +% +% tol - tolerance for stopping criterion. +% - DEFAULT 1e-7 if omitted or -1. +% maxIter - maxilambdam number of iterations +% - DEFAULT 10000, if omitted or -1. +% lineSearchFlag - 1 if line search is to be done every iteration +% - DEFAULT 0, if omitted or -1. +% continuationFlag - 1 if a continuation is to be done on the parameter lambda +% - DEFAULT 1, if omitted or -1. +% eta - line search parameter, should be in (0,1) +% - ignored if lineSearchFlag is 0. +% - DEFAULT 0.9, if omitted or -1. +% lambda - relaxation parameter +% - ignored if continuationFlag is 1. +% - DEFAULT 1e-3, if omitted or -1. +% outputFileName - Details of each iteration are dumped here, if provided. +% +% x_hat - estimate of coeeficient vector +% numIter - number of iterations until convergence +% +% +% References +% "Robust PCA: Exact Recovery of Corrupted Low-Rank Matrices via Convex Optimization", J. Wright et al., preprint 2009. +% "An Accelerated Proximal Gradient Algorithm for Nuclear Norm Regularized Least Squares problems", K.-C. Toh and S. Yun, preprint 2009. +% +% Arvind Ganesh, Summer 2009. Questions? abalasu2@illinois.edu + +DEBUG = 0 ; + +STOPPING_GROUND_TRUTH = -1; +STOPPING_DUALITY_GAP = 1; +STOPPING_SPARSE_SUPPORT = 2; +STOPPING_OBJECTIVE_VALUE = 3; +STOPPING_SUBGRADIENT = 4; +STOPPING_DEFAULT = STOPPING_SUBGRADIENT; + +stoppingCriterion = STOPPING_DEFAULT; +maxIter = 1000 ; +tolerance = 1e-3; +[m,n] = size(A) ; +x0 = zeros(n,1) ; +xG = []; + +%% Initializing optimization variables +t_k = 1 ; +t_km1 = 1 ; +L0 = 1 ; +G = A'*A ; +nIter = 0 ; +c = A'*b ; +lambda0 = 0.99*L0*norm(c,inf) ; +eta = 0.6 ; +lambda_bar = 1e-4*lambda0 ; +xk = zeros(n,1) ; +lambda = lambda0 ; +L = L0 ; +beta = 1.5 ; + +% Parse the optional inputs. +if (mod(length(varargin), 2) ~= 0 ), + error(['Extra Parameters passed to the function ''' mfilename ''' lambdast be passed in pairs.']); +end +parameterCount = length(varargin)/2; + +for parameterIndex = 1:parameterCount, + parameterName = varargin{parameterIndex*2 - 1}; + parameterValue = varargin{parameterIndex*2}; + switch lower(parameterName) + case 'stoppingcriterion' + stoppingCriterion = parameterValue; + case 'groundtruth' + xG = parameterValue; + case 'tolerance' + tolerance = parameterValue; + case 'linesearchflag' + lineSearchFlag = parameterValue; + case 'lambda' + lambda_bar = parameterValue; + case 'maxiteration' + maxIter = parameterValue; + case 'isnonnegative' + isNonnegative = parameterValue; + case 'continuationflag' + continuationFlag = parameterValue; + case 'initialization' + xk = parameterValue; + if ~all(size(xk)==[n,1]) + error('The dimension of the initial xk does not match.'); + end + case 'eta' + eta = parameterValue; + if ( eta <= 0 || eta >= 1 ) + disp('Line search parameter out of bounds, switching to default 0.9') ; + eta = 0.9 ; + end + otherwise + error(['The parameter ''' parameterName ''' is not recognized by the function ''' mfilename '''.']); + end +end +clear varargin + +if stoppingCriterion==STOPPING_GROUND_TRUTH && isempty(xG) + error('The stopping criterion must provide the ground truth value of x.'); +end + +keep_going = 1 ; +nz_x = (abs(xk)> eps*10); +f = 0.5*norm(b-A*xk)^2 + lambda_bar * norm(xk,1); +xkm1 = xk; +while keep_going && (nIter < maxIter) + nIter = nIter + 1 ; + + yk = xk + ((t_km1-1)/t_k)*(xk-xkm1) ; + + stop_backtrack = 0 ; + + temp = G*yk - c ; % gradient of f at yk + + while ~stop_backtrack + + gk = yk - (1/L)*temp ; + + xkp1 = soft(gk,lambda/L) ; + + temp1 = 0.5*norm(b-A*xkp1)^2 ; + temp2 = 0.5*norm(b-A*yk)^2 + (xkp1-yk)'*temp + (L/2)*norm(xkp1-yk)^2 ; + + if temp1 <= temp2 + stop_backtrack = 1 ; + else + L = L*beta ; + end + + end + + switch stoppingCriterion + case STOPPING_GROUND_TRUTH + keep_going = norm(xG-xkp1)>tolerance; + case STOPPING_SUBGRADIENT + sk = L*(yk-xkp1) + G*(xkp1-yk) ; + keep_going = norm(sk) > tolerance*L*max(1,norm(xkp1)); + case STOPPING_SPARSE_SUPPORT + % compute the stopping criterion based on the change + % of the number of non-zero components of the estimate + nz_x_prev = nz_x; + nz_x = (abs(xkp1)>eps*10); + num_nz_x = sum(nz_x(:)); + num_changes_active = (sum(nz_x(:)~=nz_x_prev(:))); + if num_nz_x >= 1 + criterionActiveSet = num_changes_active / num_nz_x; + keep_going = (criterionActiveSet > tolerance); + end + case STOPPING_OBJECTIVE_VALUE + % compute the stopping criterion based on the relative + % variation of the objective function. + prev_f = f; + f = 0.5*norm(b-A*xkp1)^2 + lambda_bar * norm(xk,1); + criterionObjective = abs(f-prev_f)/(prev_f); + keep_going = (criterionObjective > tolerance); + case STOPPING_DUALITY_GAP + error('Duality gap is not a valid stopping criterion for PGBP.'); + otherwise + error('Undefined stopping criterion.'); + end + + lambda = max(eta*lambda,lambda_bar) ; + + + t_kp1 = 0.5*(1+sqrt(1+4*t_k*t_k)) ; + + t_km1 = t_k ; + t_k = t_kp1 ; + xkm1 = xk ; + xk = xkp1 ; +end + +x_hat = xk ; + +function y = soft(x,T) +if sum(abs(T(:)))==0 + y = x; +else + y = max(abs(x) - T, 0); + y = sign(x).*y; +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/add_dc.m Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,33 @@ +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 + + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/addtocols.c Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,85 @@ +/************************************************************************** + * + * File name: addtocols.c + * + * Ron Rubinstein + * Computer Science Department + * Technion, Haifa 32000 Israel + * ronrubin@cs + * + * Last Updated: 19.4.2009 + * + *************************************************************************/ + + +#include "mex.h" + + +/* Input Arguments */ + +#define X_IN prhs[0] +#define V_IN prhs[1] + + +/* Output Arguments */ + +#define Y_OUT plhs[0] + + +void mexFunction(int nlhs, mxArray *plhs[], + int nrhs, const mxArray*prhs[]) + +{ + double *x, *y, *v, *xend; + mwSize m,n,m1,n1; + mwIndex counter; + + + /* Check for proper number of arguments */ + + if (nrhs != 2) { + mexErrMsgTxt("Two input arguments required."); + } else if (nlhs > 1) { + mexErrMsgTxt("Too many output arguments."); + } + + + /* Check the the input dimensions */ + + m = mxGetM(X_IN); + n = mxGetN(X_IN); + if (!mxIsDouble(X_IN) || mxIsComplex(X_IN) || mxGetNumberOfDimensions(X_IN)>2) { + mexErrMsgTxt("ADDTOCOLS requires that X be a double matrix."); + } + m1 = mxGetM(V_IN); + n1 = mxGetN(V_IN); + if (!mxIsDouble(V_IN) || mxIsComplex(V_IN) || (m1!=1 && n1!=1)) { + mexErrMsgTxt("ADDTOCOLS requires that V be a double vector."); + } + if ((m1==1 && n1!=n) || (n1==1 && m1!=n)) { + mexErrMsgTxt("Error in ADDTOCOLS: dimensions of V and X must agree."); + } + + + /* Create a matrix for the return argument */ + Y_OUT = mxCreateDoubleMatrix(m, n, mxREAL); + + + /* Assign pointers to the various parameters */ + x = mxGetPr(X_IN); + v = mxGetPr(V_IN); + y = mxGetPr(Y_OUT); + + + /* Do the actual computation */ + + xend = x+(m*n); + counter = 0; + while (x<xend) { + (*y) = (*x) + (*v); + y++; x++; counter++; + if (counter==m) {v++; counter=0;} + } + + return; +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/addtocols.m Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,13 @@ +%ADDTOCOLS Add values to the columns of a matrix. +% Y=ADDTOCOLS(X,V) adds to each column of the MxN matrix X the +% corresponding value from the N-element vector V. +% +% See also NORMCOLS. + + +% Ron Rubinstein +% Computer Science Department +% Technion, Haifa 32000 Israel +% ronrubin@cs +% +% June 2005
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/col2imstep.c Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,146 @@ +/************************************************************************** + * + * File name: col2imstep.c + * + * Ron Rubinstein + * Computer Science Department + * Technion, Haifa 32000 Israel + * ronrubin@cs + * + * Last Updated: 31.8.2009 + * + *************************************************************************/ + + +#include "mex.h" + + +/* Input Arguments */ + +#define B_IN prhs[0] +#define N_IN prhs[1] +#define SZ_IN prhs[2] +#define S_IN prhs[3] + + +/* Output Arguments */ + +#define X_OUT plhs[0] + + +void mexFunction(int nlhs, mxArray *plhs[], + int nrhs, const mxArray*prhs[]) + +{ + double *x, *b, *s; + mwSize sz[3], stepsize[3], n[3], ndims; + mwIndex i, j, k, l, m, t, blocknum; + + + /* Check for proper number of arguments */ + + if (nrhs < 3 || nrhs > 4) { + mexErrMsgTxt("Invalid number of input arguments."); + } else if (nlhs > 1) { + mexErrMsgTxt("Too many output arguments."); + } + + + /* Check the the input dimensions */ + + if (!mxIsDouble(B_IN) || mxIsComplex(B_IN) || mxGetNumberOfDimensions(B_IN)>2) { + mexErrMsgTxt("B should be a double matrix."); + } + if (!mxIsDouble(N_IN) || mxIsComplex(N_IN) || mxGetNumberOfDimensions(N_IN)>2) { + mexErrMsgTxt("Invalid output matrix size."); + } + ndims = mxGetM(N_IN)*mxGetN(N_IN); + if (ndims<2 || ndims>3) { + mexErrMsgTxt("Output matrix can only be 2-D or 3-D."); + } + if (!mxIsDouble(SZ_IN) || mxIsComplex(SZ_IN) || mxGetNumberOfDimensions(SZ_IN)>2 || mxGetM(SZ_IN)*mxGetN(SZ_IN)!=ndims) { + mexErrMsgTxt("Invalid block size."); + } + if (nrhs == 4) { + if (!mxIsDouble(S_IN) || mxIsComplex(S_IN) || mxGetNumberOfDimensions(S_IN)>2 || mxGetM(S_IN)*mxGetN(S_IN)!=ndims) { + mexErrMsgTxt("Invalid step size."); + } + } + + + /* Get parameters */ + + s = mxGetPr(N_IN); + if (s[0]<1 || s[1]<1 || (ndims==3 && s[2]<1)) { + mexErrMsgTxt("Invalid output matrix size."); + } + n[0] = (mwSize)(s[0] + 0.01); + n[1] = (mwSize)(s[1] + 0.01); + n[2] = ndims==3 ? (mwSize)(s[2] + 0.01) : 1; + + s = mxGetPr(SZ_IN); + if (s[0]<1 || s[1]<1 || (ndims==3 && s[2]<1)) { + mexErrMsgTxt("Invalid block size."); + } + sz[0] = (mwSize)(s[0] + 0.01); + sz[1] = (mwSize)(s[1] + 0.01); + sz[2] = ndims==3 ? (mwSize)(s[2] + 0.01) : 1; + + if (nrhs == 4) { + s = mxGetPr(S_IN); + if (s[0]<1 || s[1]<1 || (ndims==3 && s[2]<1)) { + mexErrMsgTxt("Invalid step size."); + } + stepsize[0] = (mwSize)(s[0] + 0.01); + stepsize[1] = (mwSize)(s[1] + 0.01); + stepsize[2] = ndims==3 ? (mwSize)(s[2] + 0.01) : 1; + } + else { + stepsize[0] = stepsize[1] = stepsize[2] = 1; + } + + if (n[0]<sz[0] || n[1]<sz[1] || (ndims==3 && n[2]<sz[2])) { + mexErrMsgTxt("Block size too large."); + } + + if (mxGetN(B_IN) != ((n[0]-sz[0])/stepsize[0]+1)*((n[1]-sz[1])/stepsize[1]+1)*((n[2]-sz[2])/stepsize[2]+1)) { + mexErrMsgTxt("Invalid number of columns in B. Please use IM2COLSTEP to compute B."); + } + + + /* Create a matrix for the return argument */ + + X_OUT = mxCreateNumericArray(ndims, n, mxDOUBLE_CLASS, mxREAL); + + + /* Assign pointers */ + + b = mxGetPr(B_IN); + x = mxGetPr(X_OUT); + + + /* Do the actual computation */ + + blocknum = 0; + + /* iterate over all blocks */ + for (k=0; k<=n[2]-sz[2]; k+=stepsize[2]) { + for (j=0; j<=n[1]-sz[1]; j+=stepsize[1]) { + for (i=0; i<=n[0]-sz[0]; i+=stepsize[0]) { + + /* add single block */ + for (m=0; m<sz[2]; m++) { + for (l=0; l<sz[1]; l++) { + for (t=0; t<sz[0]; t++) { + (x+(k+m)*n[0]*n[1]+(j+l)*n[0]+i)[t] += (b + blocknum*sz[0]*sz[1]*sz[2] + m*sz[0]*sz[1] + l*sz[0])[t]; + } + } + } + blocknum++; + + } + } + } + + return; +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/col2imstep.m Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,25 @@ +%COL2IMSTEP Rearrange matrix columns into blocks. +% A = COL2IMSTEP(B,[MM NN],[N1 N2]) rearranges the columns of B into +% sliding N1-by-N2 blocks producing the matrix A of size MM-by-NN. B is +% usually the result of calling IM2COLSTEP(...) with a stepsize of 1, or +% using Matlab's IM2COL(..,'sliding'). Overlapping blocks are summed in A. +% +% A = COL2IMSTEP(B,[MM NN],[N1 N2],[S1 S2]) arranges the blocks in A with +% a step size of (S1,S2) between them. The first block is at A(1:N1,1:N2), +% and the rest are at A((1:N1)+i*S1,(1:N2)+j*S2). Overlapping blocks are +% summed in A. Note that B is usually the result of calling +% IM2COLSTEP(...) with a stepsize of [S1 S2]. +% +% A = IM2COLSTEP(B,[MM NN KK],[N1 N2 N3],[S1 S2 S3]) generates a 3-D +% output matrix A. The step size [S1 S2 S3] may be ommitted, and defaults +% to [1 1 1]. +% +% See also IM2COLSTEP, IM2COL, COUNTCOVER. + + +% Ron Rubinstein +% Computer Science Department +% Technion, Haifa 32000 Israel +% ronrubin@cs +% +% August 2009
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/collincomb.c Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,165 @@ +/************************************************************************** + * + * File name: collincomb.c + * + * Ron Rubinstein + * Computer Science Department + * Technion, Haifa 32000 Israel + * ronrubin@cs + * + * Last Updated: 21.5.2009 + * + *************************************************************************/ + + +#include "mex.h" + + +/* Input Arguments */ + +#define A_IN prhs[0] +#define ROWS_IN prhs[1] +#define COLS_IN1 prhs[1] +#define COLS_IN2 prhs[2] +#define X_IN1 prhs[2] +#define X_IN2 prhs[3] + + +/* Output Arguments */ + +#define Y_OUT plhs[0] + + +void mexFunction(int nlhs, mxArray *plhs[], +int nrhs, const mxArray*prhs[]) + +{ + double *A, *x, *y, *rows, *cols; + mwSize m,n,m1,n1,m2,n2,rownum,colnum; + mwIndex col_id,*row_ids,i,j; + int rownumspecified=0; + + + /* Check for proper number of arguments */ + + if (nrhs!=3 && nrhs!=4) { + mexErrMsgTxt("Invalid number of arguments."); + } else if (nlhs > 1) { + mexErrMsgTxt("Too many output arguments."); + } + + + /* Check the input dimensions */ + + m = mxGetM(A_IN); + n = mxGetN(A_IN); + if (!mxIsDouble(A_IN) || mxIsComplex(A_IN) || mxGetNumberOfDimensions(A_IN)>2) { + mexErrMsgTxt("COLLINCOMB requires that A be a double matrix."); + } + + if (nrhs==3) { + + m1 = mxGetM(COLS_IN1); + n1 = mxGetN(COLS_IN1); + if (!mxIsDouble(COLS_IN1) || mxIsComplex(COLS_IN1) || (m1!=1 && n1!=1)) { + mexErrMsgTxt("COLLINCOMB requires that COLS be an index vector of type double."); + } + colnum = (m1 > n1) ? m1 : n1; /* the number of columns in the linear combination */ + + m2 = mxGetM(X_IN1); + n2 = mxGetN(X_IN1); + if (!mxIsDouble(X_IN1) || mxIsComplex(X_IN1) || (m2!=1 && n2!=1)) { + mexErrMsgTxt("COLLINCOMB requires that X be a double vector."); + } + + if (m2!=colnum && n2!=colnum) { + mexErrMsgTxt("The length of X does not match the number of columns in COLS."); + } + + rows = 0; + Y_OUT = mxCreateDoubleMatrix(m, 1, mxREAL); + cols = mxGetPr(COLS_IN1); + x = mxGetPr(X_IN1); + } + else { + + m1 = mxGetM(ROWS_IN); + n1 = mxGetN(ROWS_IN); + if (!mxIsDouble(ROWS_IN) || mxIsComplex(ROWS_IN) || (m1!=1 && n1!=1)) { + mexErrMsgTxt("COLLINCOMB requires that ROWS be an index vector of type double."); + } + rownum = (m1 > n1) ? m1 : n1; /* the number of rows in the linear combination */ + rownumspecified = 1; + rows = mxGetPr(ROWS_IN); + + m1 = mxGetM(COLS_IN2); + n1 = mxGetN(COLS_IN2); + if (!mxIsDouble(COLS_IN2) || mxIsComplex(COLS_IN2) || (m1!=1 && n1!=1)) { + mexErrMsgTxt("COLLINCOMB requires that COLS be an index vector of type double."); + } + colnum = (m1 > n1) ? m1 : n1; /* the number of columns in the linear combination */ + + m2 = mxGetM(X_IN2); + n2 = mxGetN(X_IN2); + if (!mxIsDouble(X_IN2) || mxIsComplex(X_IN2) || (m2!=1 && n2!=1)) { + mexErrMsgTxt("COLLINCOMB requires that X be a double vector."); + } + + if (m2!=colnum && n2!=colnum) { + mexErrMsgTxt("The length of X does not match the number of columns in COLS."); + } + + Y_OUT = mxCreateDoubleMatrix(rownum, 1, mxREAL); + cols = mxGetPr(COLS_IN2); + x = mxGetPr(X_IN2); + } + + + /* Assign pointers to the various parameters */ + A = mxGetPr(A_IN); + y = mxGetPr(Y_OUT); + + + if (rownumspecified) { + + /* check row indices */ + + row_ids = (mwIndex*)mxMalloc(rownum*sizeof(mwIndex)); + + for (i=0; i<rownum; ++i) { + row_ids[i] = (mwIndex)(rows[i]+0.1)-1; + if (row_ids[i]<0 || row_ids[i]>=m) { + mexErrMsgTxt("Row index in ROWS is out of range."); + } + } + + /* Do the actual computation */ + for (i=0; i<colnum; ++i) { + col_id = (mwIndex)(cols[i]+0.1)-1; + if (col_id<0 || col_id>=n) { + mexErrMsgTxt("Column index in COLS is out of range."); + } + for (j=0; j<rownum; ++j) { + y[j] += A[m*col_id+row_ids[j]]*x[i]; + } + } + + mxFree(row_ids); + } + + else { + + /* Do the actual computation */ + for (i=0; i<colnum; ++i) { + col_id = (mwIndex)(cols[i]+0.1)-1; + if (col_id<0 || col_id>=n) { + mexErrMsgTxt("Column index in COLS is out of range."); + } + for (j=0; j<m; ++j) { + y[j] += A[m*col_id+j]*x[i]; + } + } + } + + return; +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/collincomb.m Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,27 @@ +%COLLINCOMB Linear combination of matrix columns. +% Y = COLLINCOMB(A,COLS,X) computes a linear combination of the columns of +% the matrix A. The column indices are specified in the vector COLS, and +% the correspoinding coefficients are specified in the vector X. The +% vectors COLS and X must be of the same length. +% +% The call Y = COLLINCOMB(A,COLS,X) is essentially equivalent to +% +% Y = A(:,COLS)*X . +% +% However, it is implemented much more efficiently. +% +% Y = COLLINCOMB(A,ROWS,COLS,X) only works on the rows of A specified +% in ROWS, returning a vector of length equal to ROWS. This call is +% essentially equivalent to the command +% +% Y = A(ROWS,COLS)*X . +% +% See also ROWLINCOMB. + + +% Ron Rubinstein +% Computer Science Department +% Technion, Haifa 32000 Israel +% ronrubin@cs +% +% April 2009
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/completeOps.m Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,78 @@ +function data = completeOps(data) + +% Copyright 2008, Ewout van den Berg and Michael P. Friedlander +% http://www.cs.ubc.ca/labs/scl/sparco +% $Id: completeOps.m 1040 2008-06-26 20:29:02Z ewout78 $ + +operators = {}; +flagM = 0; if isfield(data,'M'), flagM = 1; end; +flagB = 0; if isfield(data,'B'), flagB = 1; end; + +if (~flagM) && (~flagB) + error('At least one of the operators M or B has be to given.'); +end + +% Define measurement matrix if needed +if ~flagM + info = data.B([],0); + data.M = opDirac(info{1}); +else + operators{end+1} = data.M; +end + +% Define sparsity basis if needed +if ~flagB + info = data.M([],0); + data.B = opDirac(info{2}); +else + operators{end+1} = data.B; +end + +% Define operator A if needed +if ~isfield(data,'A') + if (length(operators) > 1) + data.A = opFoG(operators{:}); + else + data.A = operators{1}; + end +end + +% Define empty solution if needed +if ~isfield(data,'x0') + data.x0 = []; +end + +% Define the operator size and string +opInfo = data.A([],0); +data.sizeA = [opInfo{1},opInfo{2}]; +opInfo = data.B([],0); +data.sizeB = [opInfo{1},opInfo{2}]; +opInfo = data.M([],0); +data.sizeM = [opInfo{1},opInfo{2}]; +data.op.strA = opToString(data.A); +data.op.strB = opToString(data.B); +data.op.strM = opToString(data.M); + +% Get the size of the desired signal +if ~isfield(data,'signalSize') + if ~isfield(data,'signal') + error(['At least one of the fields signal ', ... + 'or signalSize must be given.']); + end + data.signalSize = size(data.signal); +end + +% Reconstruct signal from sparse coefficients +if ~isfield(data,'reconstruct') + data.reconstruct = @(x) reshape(data.B(x,1),data.signalSize); +end + +% Reorder the fields (sort alphabetically) +m = fieldnames(data); +m = sort(m); +dataReorder = struct(); +for i=1:length(m) + eval(sprintf('dataReorder.%s = data.%s;',m{i},m{i})); +end + +data = dataReorder;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/countcover.m Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,33 @@ +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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/diag_ids.m Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,30 @@ +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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/dictdist.m Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,61 @@ +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;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/genPDF.m Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,113 @@ +function [pdf,val] = genPDF(imSize,p,pctg,distType,radius,disp) + +%[pdf,val] = genPDF(imSize,p,pctg [,distType,radius,disp]) +% +% generates a pdf for a 1d or 2d random sampling pattern +% with polynomial variable density sampling +% +% Input: +% imSize - size of matrix or vector +% p - power of polynomial +% pctg - partial sampling factor e.g. 0.5 for half +% distType - 1 or 2 for L1 or L2 distance measure +% radius - radius of fully sampled center +% disp - display output +% +% Output: +% pdf - the pdf +% val - min sampling density +% +% +% Example: +% [pdf,val] = genPDF([128,128],2,0.5,2,0,1); +% +% (c) Michael Lustig 2007 + +% This file is used with the kind permission of Michael Lustig +% (mlustig@stanford.edu), and originally appeared in the +% SparseMRI toolbox, http://www.stanford.edu/~mlustig/ . +% +% $Id: genPDF.m 1040 2008-06-26 20:29:02Z ewout78 $ + + +if nargin < 4 + distType = 2; +end + +if nargin < 5 + radius = 0; +end + +if nargin < 6 + disp = 0; +end + + +minval=0; +maxval=1; +val = 0.5; + +if length(imSize)==1 + imSize = [imSize,1]; +end + +sx = imSize(1); +sy = imSize(2); +PCTG = floor(pctg*sx*sy); + + +if sum(imSize==1)==0 % 2D + [x,y] = meshgrid(linspace(-1,1,sy),linspace(-1,1,sx)); + switch distType + case 1 + r = max(abs(x),abs(y)); + otherwise + r = sqrt(x.^2+y.^2); + r = r/max(abs(r(:))); + end + +else %1d + r = abs(linspace(-1,1,max(sx,sy))); +end + + + + +idx = find(r<radius); + +pdf = (1-r).^p; pdf(idx) = 1; +if floor(sum(pdf(:))) > PCTG + error('infeasible without undersampling dc, increase p'); +end + +% begin bisection +while(1) + val = minval/2 + maxval/2; + pdf = (1-r).^p + val; pdf(find(pdf>1)) = 1; pdf(idx)=1; + N = floor(sum(pdf(:))); + if N > PCTG % infeasible + maxval=val; + end + if N < PCTG % feasible, but not optimal + minval=val; + end + if N==PCTG % optimal + break; + end +end + +if disp + figure, + subplot(211), imshow(pdf) + if sum(imSize==1)==0 + subplot(212), plot(pdf(end/2+1,:)); + else + subplot(212), plot(pdf); + end +end + + + + + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/genSampling.m Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,53 @@ +function [minIntrVec,stat,actpctg] = genSampling(pdf,iter,tol) + +%[mask,stat,N] = genSampling(pdf,iter,tol) +% +% a monte-carlo algorithm to generate a sampling pattern with +% minimum peak interference. The number of samples will be +% sum(pdf) +- tol +% +% pdf - probability density function to choose samples from +% iter - number of tries +% tol - the deviation from the desired number of samples in samples +% +% returns: +% mask - sampling pattern +% stat - vector of min interferences measured each try +% actpctg - actual undersampling factor +% +% (c) Michael Lustig 2007 + +% This file is used with the kind permission of Michael Lustig +% (mlustig@stanford.edu), and originally appeared in the +% SparseMRI toolbox, http://www.stanford.edu/~mlustig/ . +% +% $Id: genSampling.m 1040 2008-06-26 20:29:02Z ewout78 $ + +% h = waitbar(0); + +pdf(find(pdf>1)) = 1; +K = sum(pdf(:)); + +minIntr = 1e99; +minIntrVec = zeros(size(pdf)); + +for n=1:iter + tmp = zeros(size(pdf)); + while abs(sum(tmp(:)) - K) > tol + tmp = rand(size(pdf))<pdf; + end + + TMP = ifft2(tmp./pdf); + if max(abs(TMP(2:end))) < minIntr + minIntr = max(abs(TMP(2:end))); + minIntrVec = tmp; + end + stat(n) = max(abs(TMP(2:end))); + % waitbar(n/iter,h); +end + +actpctg = sum(minIntrVec(:))/prod(size(minIntrVec)); + +% close(h); + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/im2colstep.c Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,131 @@ +/************************************************************************** + * + * File name: im2colstep.c + * + * Ron Rubinstein + * Computer Science Department + * Technion, Haifa 32000 Israel + * ronrubin@cs + * + * Last Updated: 31.8.2009 + * + *************************************************************************/ + + +#include "mex.h" +#include <string.h> + + +/* Input Arguments */ + +#define X_IN prhs[0] +#define SZ_IN prhs[1] +#define S_IN prhs[2] + + +/* Output Arguments */ + +#define B_OUT plhs[0] + + +void mexFunction(int nlhs, mxArray *plhs[], + int nrhs, const mxArray*prhs[]) + +{ + double *x, *b, *s; + mwSize sz[3], stepsize[3], n[3], ndims; + mwIndex i, j, k, l, m, blocknum; + + + /* Check for proper number of arguments */ + + if (nrhs < 2 || nrhs > 3) { + mexErrMsgTxt("Invalid number of input arguments."); + } else if (nlhs > 1) { + mexErrMsgTxt("Too many output arguments."); + } + + + /* Check the the input dimensions */ + + ndims = mxGetNumberOfDimensions(X_IN); + + if (!mxIsDouble(X_IN) || mxIsComplex(X_IN) || ndims>3) { + mexErrMsgTxt("X should be a 2-D or 3-D double matrix."); + } + if (!mxIsDouble(SZ_IN) || mxIsComplex(SZ_IN) || mxGetNumberOfDimensions(SZ_IN)>2 || mxGetM(SZ_IN)*mxGetN(SZ_IN)!=ndims) { + mexErrMsgTxt("Invalid block size."); + } + if (nrhs == 3) { + if (!mxIsDouble(S_IN) || mxIsComplex(S_IN) || mxGetNumberOfDimensions(S_IN)>2 || mxGetM(S_IN)*mxGetN(S_IN)!=ndims) { + mexErrMsgTxt("Invalid step size."); + } + } + + + /* Get parameters */ + + s = mxGetPr(SZ_IN); + if (s[0]<1 || s[1]<1 || (ndims==3 && s[2]<1)) { + mexErrMsgTxt("Invalid block size."); + } + sz[0] = (mwSize)(s[0] + 0.01); + sz[1] = (mwSize)(s[1] + 0.01); + sz[2] = ndims==3 ? (mwSize)(s[2] + 0.01) : 1; + + if (nrhs == 3) { + s = mxGetPr(S_IN); + if (s[0]<1 || s[1]<1 || (ndims==3 && s[2]<1)) { + mexErrMsgTxt("Invalid step size."); + } + stepsize[0] = (mwSize)(s[0] + 0.01); + stepsize[1] = (mwSize)(s[1] + 0.01); + stepsize[2] = ndims==3 ? (mwSize)(s[2] + 0.01) : 1; + } + else { + stepsize[0] = stepsize[1] = stepsize[2] = 1; + } + + n[0] = (mxGetDimensions(X_IN))[0]; + n[1] = (mxGetDimensions(X_IN))[1]; + n[2] = ndims==3 ? (mxGetDimensions(X_IN))[2] : 1; + + if (n[0]<sz[0] || n[1]<sz[1] || (ndims==3 && n[2]<sz[2])) { + mexErrMsgTxt("Block size too large."); + } + + + /* Create a matrix for the return argument */ + + B_OUT = mxCreateDoubleMatrix(sz[0]*sz[1]*sz[2], ((n[0]-sz[0])/stepsize[0]+1)*((n[1]-sz[1])/stepsize[1]+1)*((n[2]-sz[2])/stepsize[2]+1), mxREAL); + + + /* Assign pointers */ + + x = mxGetPr(X_IN); + b = mxGetPr(B_OUT); + + + /* Do the actual computation */ + + blocknum = 0; + + /* iterate over all blocks */ + for (k=0; k<=n[2]-sz[2]; k+=stepsize[2]) { + for (j=0; j<=n[1]-sz[1]; j+=stepsize[1]) { + for (i=0; i<=n[0]-sz[0]; i+=stepsize[0]) { + + /* copy single block */ + for (m=0; m<sz[2]; m++) { + for (l=0; l<sz[1]; l++) { + memcpy(b + blocknum*sz[0]*sz[1]*sz[2] + m*sz[0]*sz[1] + l*sz[0], x+(k+m)*n[0]*n[1]+(j+l)*n[0]+i, sz[0]*sizeof(double)); + } + } + blocknum++; + + } + } + } + + return; +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/im2colstep.m Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,31 @@ +%IM2COLSTEP Rearrange matrix blocks into columns. +% B = IM2COLSTEP(A,[N1 N2]) converts each sliding N1-by-N2 block of the +% 2-D matrix A into a column of B, with no zero padding. B has N1*N2 rows +% and will contain as many columns as there are N1-by-N2 neighborhoods in +% A. Each column of B contains a neighborhood of A reshaped as NHOOD(:), +% where NHOOD is a matrix containing an N1-by-N2 neighborhood of A. +% +% B = IM2COLSTEP(A,[N1 N2],[S1 S2]) extracts neighborhoods of A with a +% step size of (S1,S2) between them. The first extracted neighborhood is +% A(1:N1,1:N2), and the rest are of the form A((1:N1)+i*S1,(1:N2)+j*S2). +% Note that to ensure coverage of all A by neighborhoods, +% (size(A,i)-Ni)/Si must be whole for i=1,2. The default function behavior +% corresponds to [S1 S2] = [1 1]. Setting S1>=N1 and S2>=N2 results in no +% overlap between the neighborhoods. +% +% B = IM2COLSTEP(A,[N1 N2 N3],[S1 S2 S3]) operates on a 3-D matrix A. The +% step size [S1 S2 S3] may be ommitted, and defaults to [1 1 1]. +% +% Note: the call IM2COLSTEP(A,[N1 N2]) produces the same output as +% Matlab's IM2COL(A,[N1 N2],'sliding'). However, it is significantly +% faster. +% +% See also COL2IMSTEP, IM2COL, COUNTCOVER. + + +% Ron Rubinstein +% Computer Science Department +% Technion, Haifa 32000 Israel +% ronrubin@cs +% +% August 2009
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/imnormalize.m Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,19 @@ +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);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/iswhole.m Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,17 @@ +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;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/make.m Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,40 @@ +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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/mexutils.c Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,79 @@ +/************************************************************************** + * + * File name: mexutils.c + * + * Ron Rubinstein + * Computer Science Department + * Technion, Haifa 32000 Israel + * ronrubin@cs + * + * Last Updated: 15.8.2009 + * + *************************************************************************/ + +#include "mexutils.h" +#include <math.h> + + + +/* verify that the mxArray contains a double matrix */ + +void checkmatrix(const mxArray *param, char *fname, char *pname) +{ + char errmsg[100]; + sprintf(errmsg, "%.15s requires that %.25s be a double matrix.", fname, pname); + if (!mxIsDouble(param) || mxIsComplex(param) || mxGetNumberOfDimensions(param)>2) { + mexErrMsgTxt(errmsg); + } +} + + +/* verify that the mxArray contains a 1-D double vector */ + +void checkvector(const mxArray *param, char *fname, char *pname) +{ + char errmsg[100]; + sprintf(errmsg, "%.15s requires that %.25s be a double vector.", fname, pname); + if (!mxIsDouble(param) || mxIsComplex(param) || mxGetNumberOfDimensions(param)>2 || (mxGetM(param)!=1 && mxGetN(param)!=1)) { + mexErrMsgTxt(errmsg); + } +} + + +/* verify that the mxArray contains a double scalar */ + +void checkscalar(const mxArray *param, char *fname, char *pname) +{ + char errmsg[100]; + sprintf(errmsg, "%.15s requires that %.25s be a double scalar.", fname, pname); + if (!mxIsDouble(param) || mxIsComplex(param) || mxGetNumberOfDimensions(param)>2 || + mxGetM(param)!=1 || mxGetN(param)!=1) + { + mexErrMsgTxt(errmsg); + } +} + + +/* verify that the mxArray contains a sparse matrix */ + +void checksparse(const mxArray *param, char *fname, char *pname) +{ + char errmsg[100]; + sprintf(errmsg, "%.15s requires that %.25s be sparse.", fname, pname); + if (!mxIsSparse(param)) { + mexErrMsgTxt(errmsg); + } +} + + +/* verify that the mxArray contains a 1-dimensional cell array */ + +void checkcell_1d(const mxArray *param, char *fname, char *pname) +{ + char errmsg[100]; + sprintf(errmsg, "%.15s requires that %.25s be a 1-D cell array.", fname, pname); + if (!mxIsCell(param) || mxGetNumberOfDimensions(param)>2 || (mxGetM(param)!=1 && mxGetN(param)!=1)) { + mexErrMsgTxt(errmsg); + } +} +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/mexutils.h Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,103 @@ +/************************************************************************** + * + * File name: mexutils.h + * + * Ron Rubinstein + * Computer Science Department + * Technion, Haifa 32000 Israel + * ronrubin@cs + * + * Last Updated: 18.8.2009 + * + * Utility functions for MEX files. + * + *************************************************************************/ + + +#ifndef __MEX_UTILS_H__ +#define __MEX_UTILS_H__ + +#include "mex.h" + + + +/************************************************************************** + * Function checkmatrix: + * + * Verify that the specified mxArray is real, of type double, and has + * no more than two dimensions. If not, an error message is printed + * and the mex file terminates. + * + * Parameters: + * param - the mxArray to be checked + * fname - the name of the function where the error occured (15 characters or less) + * pname - the name of the parameter (25 characters or less) + * + **************************************************************************/ +void checkmatrix(const mxArray *param, char *fname, char *pname); + + +/************************************************************************** + * Function checkvector: + * + * Verify that the specified mxArray is 1-D, real, and of type double. The + * vector may be a column or row vector. Otherwise, an error message is + * printed and the mex file terminates. + * + * Parameters: + * param - the mxArray to be checked + * fname - the name of the function where the error occured (15 characters or less) + * pname - the name of the parameter (25 characters or less) + * + **************************************************************************/ +void checkvector(const mxArray *param, char *fname, char *pname); + + +/************************************************************************** + * Function checkscalar: + * + * Verify that the specified mxArray represents a real double scalar value. + * If not, an error message is printed and the mex file terminates. + * + * Parameters: + * param - the mxArray to be checked + * fname - the name of the function where the error occured (15 characters or less) + * pname - the name of the parameter (25 characters or less) + * + **************************************************************************/ +void checkscalar(const mxArray *param, char *fname, char *pname); + + +/************************************************************************** + * Function checksparse: + * + * Verify that the specified mxArray contains a sparse matrix. If not, + * an error message is printed and the mex file terminates. + * + * Parameters: + * param - the mxArray to be checked + * fname - the name of the function where the error occured (15 characters or less) + * pname - the name of the parameter (25 characters or less) + * + **************************************************************************/ +void checksparse(const mxArray *param, char *fname, char *pname); + + +/************************************************************************** + * Function checkcell_1d: + * + * Verify that the specified mxArray is a 1-D cell array. The cell array + * may be arranged as either a column or a row. If not, an error message + * is printed and the mex file terminates. + * + * Parameters: + * param - the mxArray to be checked + * fname - the name of the function where the error occured (15 characters or less) + * pname - the name of the parameter (25 characters or less) + * + **************************************************************************/ +void checkcell_1d(const mxArray *param, char *fname, char *pname); + + +#endif +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/myblas.c Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,663 @@ +/************************************************************************** + * + * File name: myblas.c + * + * Ron Rubinstein + * Computer Science Department + * Technion, Haifa 32000 Israel + * ronrubin@cs + * + * Version: 1.1 + * Last updated: 13.8.2009 + * + *************************************************************************/ + + +#include "myblas.h" +#include <ctype.h> + + +/* find maximum of absolute values */ + +mwIndex maxabs(double c[], mwSize m) +{ + mwIndex maxid=0, k; + double absval, maxval = SQR(*c); /* use square which is quicker than absolute value */ + + for (k=1; k<m; ++k) { + absval = SQR(c[k]); + if (absval > maxval) { + maxval = absval; + maxid = k; + } + } + return maxid; +} + + +/* compute y := alpha*x + y */ + +void vec_sum(double alpha, double x[], double y[], mwSize n) +{ + mwIndex i; + + for (i=0; i<n; ++i) { + y[i] += alpha*x[i]; + } +} + + +/* compute y := alpha*A*x */ + +void mat_vec(double alpha, double A[], double x[], double y[], mwSize n, mwSize m) +{ + mwIndex i, j, i_n; + double *Ax; + + Ax = mxCalloc(n,sizeof(double)); + + for (i=0; i<m; ++i) { + i_n = i*n; + for (j=0; j<n; ++j) { + Ax[j] += A[i_n+j] * x[i]; + } + } + + for (j=0; j<n; ++j) { + y[j] = alpha*Ax[j]; + } + + mxFree(Ax); +} + + +/* compute y := alpha*A'*x */ + +void matT_vec(double alpha, double A[], double x[], double y[], mwSize n, mwSize m) +{ + mwIndex i, j, n_i; + double sum0, sum1, sum2, sum3; + + for (j=0; j<m; ++j) { + y[j] = 0; + } + + /* use loop unrolling to accelerate computation */ + + for (i=0; i<m; ++i) { + n_i = n*i; + sum0 = sum1 = sum2 = sum3 = 0; + for (j=0; j+4<n; j+=4) { + sum0 += A[n_i+j]*x[j]; + sum1 += A[n_i+j+1]*x[j+1]; + sum2 += A[n_i+j+2]*x[j+2]; + sum3 += A[n_i+j+3]*x[j+3]; + } + y[i] += alpha * ((sum0 + sum1) + (sum2 + sum3)); + while (j<n) { + y[i] += alpha*A[n_i+j]*x[j]; + j++; + } + } +} + + +/* compute y := alpha*A*x */ + +void mat_sp_vec(double alpha, double pr[], mwIndex ir[], mwIndex jc[], double x[], double y[], mwSize n, mwSize m) +{ + + mwIndex i, j, j1, j2; + + for (i=0; i<n; ++i) { + y[i] = 0; + } + + j2 = jc[0]; + for (i=0; i<m; ++i) { + j1 = j2; j2 = jc[i+1]; + for (j=j1; j<j2; ++j) { + y[ir[j]] += alpha * pr[j] * x[i]; + } + } + +} + + +/* compute y := alpha*A'*x */ + +void matT_sp_vec(double alpha, double pr[], mwIndex ir[], mwIndex jc[], double x[], double y[], mwSize n, mwSize m) +{ + + mwIndex i, j, j1, j2; + + for (i=0; i<m; ++i) { + y[i] = 0; + } + + j2 = jc[0]; + for (i=0; i<m; ++i) { + j1 = j2; j2 = jc[i+1]; + for (j=j1; j<j2; ++j) { + y[i] += alpha * pr[j] * x[ir[j]]; + } + } + +} + + +/* compute y := alpha*A*x */ + +void mat_vec_sp(double alpha, double A[], double pr[], mwIndex ir[], mwIndex jc[], double y[], mwSize n, mwSize m) +{ + + mwIndex i, j, j_n, k, kend; + + for (i=0; i<n; ++i) { + y[i] = 0; + } + + kend = jc[1]; + if (kend==0) { /* x is empty */ + return; + } + + for (k=0; k<kend; ++k) { + j = ir[k]; + j_n = j*n; + for (i=0; i<n; ++i) { + y[i] += alpha * A[i+j_n] * pr[k]; + } + } + +} + + +/* compute y := alpha*A'*x */ + +void matT_vec_sp(double alpha, double A[], double pr[], mwIndex ir[], mwIndex jc[], double y[], mwSize n, mwSize m) +{ + + mwIndex i, j, j_n, k, kend; + + for (i=0; i<m; ++i) { + y[i] = 0; + } + + kend = jc[1]; + if (kend==0) { /* x is empty */ + return; + } + + for (j=0; j<m; ++j) { + j_n = j*n; + for (k=0; k<kend; ++k) { + i = ir[k]; + y[j] += alpha * A[i+j_n] * pr[k]; + } + } + +} + + +/* compute y := alpha*A*x */ + +void mat_sp_vec_sp(double alpha, double pr[], mwIndex ir[], mwIndex jc[], double prx[], mwIndex irx[], mwIndex jcx[], double y[], mwSize n, mwSize m) +{ + + mwIndex i, j, k, kend, j1, j2; + + for (i=0; i<n; ++i) { + y[i] = 0; + } + + kend = jcx[1]; + if (kend==0) { /* x is empty */ + return; + } + + for (k=0; k<kend; ++k) { + i = irx[k]; + j1 = jc[i]; j2 = jc[i+1]; + for (j=j1; j<j2; ++j) { + y[ir[j]] += alpha * pr[j] * prx[k]; + } + } + +} + + +/* compute y := alpha*A'*x */ + +void matT_sp_vec_sp(double alpha, double pr[], mwIndex ir[], mwIndex jc[], double prx[], mwIndex irx[], mwIndex jcx[], double y[], mwSize n, mwSize m) +{ + + mwIndex i, j, k, jend, kend, jadd, kadd, delta; + + for (i=0; i<m; ++i) { + y[i] = 0; + } + + kend = jcx[1]; + if (kend==0) { /* x is empty */ + return; + } + + for (i=0; i<m; ++i) { + j = jc[i]; + jend = jc[i+1]; + k = 0; + while (j<jend && k<kend) { + + delta = ir[j] - irx[k]; + + if (delta) { /* if indices differ - increment the smaller one */ + jadd = delta<0; + kadd = 1-jadd; + j += jadd; + k += kadd; + } + + else { /* indices are equal - add to result and increment both */ + y[i] += alpha * pr[j] * prx[k]; + j++; k++; + } + } + } + +} + + +/* matrix-matrix multiplication */ + +void mat_mat(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k) +{ + mwIndex i1, i2, i3, iX, iA, i2_n; + double b; + + for (i1=0; i1<n*k; i1++) { + X[i1] = 0; + } + + for (i2=0; i2<m; ++i2) { + i2_n = i2*n; + iX = 0; + for (i3=0; i3<k; ++i3) { + iA = i2_n; + b = B[i2+i3*m]; + for (i1=0; i1<n; ++i1) { + X[iX++] += A[iA++]*b; + } + } + } + + for (i1=0; i1<n*k; i1++) { + X[i1] *= alpha; + } +} + + +/* matrix-transpose-matrix multiplication */ + +void matT_mat(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k) +{ + mwIndex i1, i2, i3, iX, iA, i2_n; + double *x, sum0, sum1, sum2, sum3; + + for (i2=0; i2<m; ++i2) { + for (i3=0; i3<k; ++i3) { + sum0 = sum1 = sum2 = sum3 = 0; + for (i1=0; i1+4<n; i1+=4) { + sum0 += A[i1+0+i2*n]*B[i1+0+i3*n]; + sum1 += A[i1+1+i2*n]*B[i1+1+i3*n]; + sum2 += A[i1+2+i2*n]*B[i1+2+i3*n]; + sum3 += A[i1+3+i2*n]*B[i1+3+i3*n]; + } + X[i2+i3*m] = (sum0+sum1) + (sum2+sum3); + while(i1<n) { + X[i2+i3*m] += A[i1+i2*n]*B[i1+i3*n]; + i1++; + } + } + } + + for (i1=0; i1<m*k; i1++) { + X[i1] *= alpha; + } +} + + +/* tensor-matrix product */ + +void tens_mat(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k, mwSize l) +{ + mwIndex i1, i2, i3, i4, i2_n, nml; + double b; + + nml = n*m*l; + for (i1=0; i1<nml; ++i1) { + X[i1] = 0; + } + + for (i2=0; i2<m; ++i2) { + i2_n = i2*n; + for (i3=0; i3<k; ++i3) { + for (i4=0; i4<l; ++i4) { + b = B[i4+i3*l]; + for (i1=0; i1<n; ++i1) { + X[i1 + i2_n + i4*n*m] += A[i1 + i2_n + i3*n*m] * b; + } + } + } + } + + for (i1=0; i1<nml; ++i1) { + X[i1] *= alpha; + } +} + + +/* tensor-matrix-transpose product */ + +void tens_matT(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k, mwSize l) +{ + mwIndex i1, i2, i3, i4, i2_n, nml; + double b; + + nml = n*m*l; + for (i1=0; i1<nml; ++i1) { + X[i1] = 0; + } + + for (i2=0; i2<m; ++i2) { + i2_n = i2*n; + for (i4=0; i4<l; ++i4) { + for (i3=0; i3<k; ++i3) { + b = B[i3+i4*k]; + for (i1=0; i1<n; ++i1) { + X[i1 + i2_n + i4*n*m] += A[i1 + i2_n + i3*n*m] * b; + } + } + } + } + + for (i1=0; i1<nml; ++i1) { + X[i1] *= alpha; + } +} + + +/* dot product */ + +double dotprod(double a[], double b[], mwSize n) +{ + double sum = 0; + mwIndex i; + for (i=0; i<n; ++i) + sum += a[i]*b[i]; + return sum; +} + + +/* find maximum of vector */ + +mwIndex maxpos(double c[], mwSize m) +{ + mwIndex maxid=0, k; + double val, maxval = *c; + + for (k=1; k<m; ++k) { + val = c[k]; + if (val > maxval) { + maxval = val; + maxid = k; + } + } + return maxid; +} + + +/* solve L*x = b */ + +void backsubst_L(double L[], double b[], double x[], mwSize n, mwSize k) +{ + mwIndex i, j; + double rhs; + + for (i=0; i<k; ++i) { + rhs = b[i]; + for (j=0; j<i; ++j) { + rhs -= L[j*n+i]*x[j]; + } + x[i] = rhs/L[i*n+i]; + } +} + + +/* solve L'*x = b */ + +void backsubst_Lt(double L[], double b[], double x[], mwSize n, mwSize k) +{ + mwIndex i, j; + double rhs; + + for (i=k; i>=1; --i) { + rhs = b[i-1]; + for (j=i; j<k; ++j) { + rhs -= L[(i-1)*n+j]*x[j]; + } + x[i-1] = rhs/L[(i-1)*n+i-1]; + } +} + + +/* solve U*x = b */ + +void backsubst_U(double U[], double b[], double x[], mwSize n, mwSize k) +{ + mwIndex i, j; + double rhs; + + for (i=k; i>=1; --i) { + rhs = b[i-1]; + for (j=i; j<k; ++j) { + rhs -= U[j*n+i-1]*x[j]; + } + x[i-1] = rhs/U[(i-1)*n+i-1]; + } +} + + +/* solve U'*x = b */ + +void backsubst_Ut(double U[], double b[], double x[], mwSize n, mwSize k) +{ + mwIndex i, j; + double rhs; + + for (i=0; i<k; ++i) { + rhs = b[i]; + for (j=0; j<i; ++j) { + rhs -= U[i*n+j]*x[j]; + } + x[i] = rhs/U[i*n+i]; + } +} + + +/* back substitution solver */ + +void backsubst(char ul, double A[], double b[], double x[], mwSize n, mwSize k) +{ + if (tolower(ul) == 'u') { + backsubst_U(A, b, x, n, k); + } + else if (tolower(ul) == 'l') { + backsubst_L(A, b, x, n, k); + } + else { + mexErrMsgTxt("Invalid triangular matrix type: must be ''U'' or ''L''"); + } +} + + +/* solve equation set using cholesky decomposition */ + +void cholsolve(char ul, double A[], double b[], double x[], mwSize n, mwSize k) +{ + double *tmp; + + tmp = mxMalloc(k*sizeof(double)); + + if (tolower(ul) == 'l') { + backsubst_L(A, b, tmp, n, k); + backsubst_Lt(A, tmp, x, n, k); + } + else if (tolower(ul) == 'u') { + backsubst_Ut(A, b, tmp, n, k); + backsubst_U(A, tmp, x, n, k); + } + else { + mexErrMsgTxt("Invalid triangular matrix type: must be either ''U'' or ''L''"); + } + + mxFree(tmp); +} + + +/* perform a permutation assignment y := x(ind(1:k)) */ + +void vec_assign(double y[], double x[], mwIndex ind[], mwSize k) +{ + mwIndex i; + + for (i=0; i<k; ++i) + y[i] = x[ind[i]]; +} + + +/* matrix transpose */ + +void transpose(double X[], double Y[], mwSize n, mwSize m) +{ + mwIndex i, j, i_m, j_n; + + if (n<m) { + for (j=0; j<m; ++j) { + j_n = j*n; + for (i=0; i<n; ++i) { + Y[j+i*m] = X[i+j_n]; + } + } + } + else { + for (i=0; i<n; ++i) { + i_m = i*m; + for (j=0; j<m; ++j) { + Y[j+i_m] = X[i+j*n]; + } + } + } +} + + +/* print contents of matrix */ + +void printmat(double A[], int n, int m, char* matname) +{ + int i, j; + mexPrintf("\n%s = \n\n", matname); + + if (n*m==0) { + mexPrintf(" Empty matrix: %d-by-%d\n\n", n, m); + return; + } + + for (i=0; i<n; ++i) { + for (j=0; j<m; ++j) + mexPrintf(" %lf", A[j*n+i]); + mexPrintf("\n"); + } + mexPrintf("\n"); +} + + +/* print contents of sparse matrix */ + +void printspmat(mxArray *a, char* matname) +{ + mwIndex *aJc = mxGetJc(a); + mwIndex *aIr = mxGetIr(a); + double *aPr = mxGetPr(a); + + int i; + + mexPrintf("\n%s = \n\n", matname); + + for (i=0; i<aJc[1]; ++i) + printf(" (%d,1) = %lf\n", aIr[i]+1,aPr[i]); + + mexPrintf("\n"); +} + + + +/* matrix multiplication using Winograd's algorithm */ + +/* +void mat_mat2(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k) +{ + + mwIndex i1, i2, i3, iX, iA, i2_n; + double b, *AA, *BB; + + AA = mxCalloc(n,sizeof(double)); + BB = mxCalloc(k,sizeof(double)); + + for (i1=0; i1<n*k; i1++) { + X[i1] = 0; + } + + for (i1=0; i1<n; ++i1) { + for (i2=0; i2<m/2; ++i2) { + AA[i1] += A[i1+2*i2*n]*A[i1+(2*i2+1)*n]; + } + } + + for (i2=0; i2<k; ++i2) { + for (i1=0; i1<m/2; ++i1) { + BB[i2] += B[2*i1+i2*m]*B[2*i1+1+i2*m]; + } + } + + for (i2=0; i2<k; ++i2) { + for (i3=0; i3<m/2; ++i3) { + for (i1=0; i1<n; ++i1) { + X[i1+i2*n] += (A[i1+(2*i3)*n]+B[2*i3+1+i2*m])*(A[i1+(2*i3+1)*n]+B[2*i3+i2*m]); + } + } + } + + if (m%2) { + for (i2=0; i2<k; ++i2) { + for (i1=0; i1<n; ++i1) { + X[i1+i2*n] += A[i1+(m-1)*n]*B[m-1+i2*m]; + } + } + } + + for (i2=0; i2<k; ++i2) { + for (i1=0; i1<n; ++i1) { + X[i1+i2*n] -= (AA[i1] + BB[i2]); + X[i1+i2*n] *= alpha; + } + } + + mxFree(AA); + mxFree(BB); +} +*/ + + + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/myblas.h Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,492 @@ +/************************************************************************** + * + * File name: myblas.h + * + * Ron Rubinstein + * Computer Science Department + * Technion, Haifa 32000 Israel + * ronrubin@cs + * + * Version: 1.1 + * Last updated: 17.8.2009 + * + * A collection of basic linear algebra functions, in the spirit of the + * BLAS/LAPACK libraries. + * + *************************************************************************/ + + + +#ifndef __MY_BLAS_H__ +#define __MY_BLAS_H__ + + +#include "mex.h" +#include <math.h> + + + +/************************************************************************** + * Squared value. + **************************************************************************/ +#define SQR(X) ((X)*(X)) + + + +/************************************************************************** + * Matrix-vector multiplication. + * + * Computes an operation of the form: + * + * y := alpha*A*x + * + * Parameters: + * A - matrix of size n X m + * x - vector of length m + * y - output vector of length n + * alpha - real constant + * n, m - dimensions of A + * + * Note: This function re-writes the contents of y. + * + **************************************************************************/ +void mat_vec(double alpha, double A[], double x[], double y[], mwSize n, mwSize m); + + + +/************************************************************************** + * Matrix-transpose-vector multiplication. + * + * Computes an operation of the form: + * + * y := alpha*A'*x + * + * Parameters: + * A - matrix of size n X m + * x - vector of length n + * y - output vector of length m + * alpha - real constant + * n, m - dimensions of A + * + * Note: This function re-writes the contents of y. + * + **************************************************************************/ +void matT_vec(double alpha, double A[], double x[], double y[], mwSize n, mwSize m); + + + +/************************************************************************** + * Sparse-matrix-vector multiplication. + * + * Computes an operation of the form: + * + * y := alpha*A*x + * + * where A is a sparse matrix. + * + * Parameters: + * pr,ir,jc - sparse representation of the matrix A, of size n x m + * x - vector of length m + * y - output vector of length n + * alpha - real constant + * n, m - dimensions of A + * + * Note: This function re-writes the contents of y. + * + **************************************************************************/ +void mat_sp_vec(double alpha, double pr[], mwIndex ir[], mwIndex jc[], double x[], double y[], mwSize n, mwSize m); + + + +/************************************************************************** + * Sparse-matrix-transpose-vector multiplication. + * + * Computes an operation of the form: + * + * y := alpha*A'*x + * + * where A is a sparse matrix. + * + * Parameters: + * pr,ir,jc - sparse representation of the matrix A, of size n x m + * x - vector of length m + * y - output vector of length n + * alpha - real constant + * n, m - dimensions of A + * + * Note: This function re-writes the contents of y. + * + **************************************************************************/ +void matT_sp_vec(double alpha, double pr[], mwIndex ir[], mwIndex jc[], double x[], double y[], mwSize n, mwSize m); + + + +/************************************************************************** + * Matrix-sparse-vector multiplication. + * + * Computes an operation of the form: + * + * y := alpha*A*x + * + * where A is a matrix and x is a sparse vector. + * + * Parameters: + * A - matrix of size n X m + * pr,ir,jc - sparse representation of the vector x, of length m + * y - output vector of length n + * alpha - real constant + * n, m - dimensions of A + * + * Note: This function re-writes the contents of y. + * + **************************************************************************/ +void mat_vec_sp(double alpha, double A[], double pr[], mwIndex ir[], mwIndex jc[], double y[], mwSize n, mwSize m); + + + +/************************************************************************** + * Matrix-transpose-sparse-vector multiplication. + * + * Computes an operation of the form: + * + * y := alpha*A'*x + * + * where A is a matrix and x is a sparse vector. + * + * Parameters: + * A - matrix of size n X m + * pr,ir,jc - sparse representation of the vector x, of length n + * y - output vector of length m + * alpha - real constant + * n, m - dimensions of A + * + * Note: This function re-writes the contents of y. + * + **************************************************************************/ +void matT_vec_sp(double alpha, double A[], double pr[], mwIndex ir[], mwIndex jc[], double y[], mwSize n, mwSize m); + + + +/************************************************************************** + * Sparse-matrix-sparse-vector multiplication. + * + * Computes an operation of the form: + * + * y := alpha*A*x + * + * where A is a sparse matrix and x is a sparse vector. + * + * Parameters: + * pr,ir,jc - sparse representation of the matrix A, of size n x m + * prx,irx,jcx - sparse representation of the vector x (of length m) + * y - output vector of length n + * alpha - real constant + * n, m - dimensions of A + * + * Note: This function re-writes the contents of y. + * + **************************************************************************/ +void mat_sp_vec_sp(double alpha, double pr[], mwIndex ir[], mwIndex jc[], double prx[], mwIndex irx[], mwIndex jcx[], double y[], mwSize n, mwSize m); + + + +/************************************************************************** + * Sparse-matrix-transpose-sparse-vector multiplication. + * + * Computes an operation of the form: + * + * y := alpha*A'*x + * + * where A is a sparse matrix and x is a sparse vector. + * + * Importnant note: this function is provided for completeness, but is NOT efficient. + * If possible, convert x to non-sparse representation and use matT_vec_sp instead. + * + * Parameters: + * pr,ir,jc - sparse representation of the matrix A, of size n x m + * prx,irx,jcx - sparse representation of the vector x (of length n) + * y - output vector of length n + * alpha - real constant + * n, m - dimensions of A + * + * Note: This function re-writes the contents of y. + * + **************************************************************************/ +void matT_sp_vec_sp(double alpha, double pr[], mwIndex ir[], mwIndex jc[], double prx[], mwIndex irx[], mwIndex jcx[], double y[], mwSize n, mwSize m); + + + +/************************************************************************** + * Matrix-matrix multiplication. + * + * Computes an operation of the form: + * + * X := alpha*A*B + * + * Parameters: + * A - matrix of size n X m + * B - matrix of size m X k + * X - output matrix of size n X k + * alpha - real constant + * n, m, k - dimensions of A, B + * + * Note: This function re-writes the contents of X. + * + **************************************************************************/ +void mat_mat(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k); + + + +/************************************************************************** + * Matrix-transpose-matrix multiplication. + * + * Computes an operation of the form: + * + * X := alpha*A*B + * + * Parameters: + * A - matrix of size n X m + * B - matrix of size m X k + * X - output matrix of size n X k + * alpha - real constant + * n, m, k - dimensions of A, B + * + * Note: This function re-writes the contents of X. + * + **************************************************************************/ +void matT_mat(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k); + + + +/************************************************************************** + * Tensor-matrix multiplication. + * + * This function accepts a 3-D tensor A of size n X m X k + * and a 2-D matrix B of size l X k. + * The function computes the 3-D tensor X of size n X m X l, where + * + * X(i,j,:) = B*A(i,j,:) + * + * for all i,j. + * + * Parameters: + * A - tensor of size n X m X k + * B - matrix of size l X k + * X - output tensor of size n X m X l + * alpha - real constant + * n, m, k, l - dimensions of A, B + * + * Note: This function re-writes the contents of X. + * + **************************************************************************/ +void tens_mat(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k, mwSize l); + + + +/************************************************************************** + * Tensor-matrix-transpose multiplication. + * + * This function accepts a 3-D tensor A of size n X m X k + * and a 2-D matrix B of size k X l. + * The function computes the 3-D tensor X of size n X m X l, where + * + * X(i,j,:) = B'*A(i,j,:) + * + * for all i,j. + * + * Parameters: + * A - tensor of size n X m X k + * B - matrix of size k X l + * X - output tensor of size n X m X l + * alpha - real constant + * n, m, k, l - dimensions of A, B + * + * Note: This function re-writes the contents of X. + * + **************************************************************************/ +void tens_matT(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k, mwSize l); + + + +/************************************************************************** + * Vector-vector sum. + * + * Computes an operation of the form: + * + * y := alpha*x + y + * + * Parameters: + * x - vector of length n + * y - output vector of length n + * alpha - real constant + * n - length of x,y + * + * Note: This function re-writes the contents of y. + * + **************************************************************************/ +void vec_sum(double alpha, double x[], double y[], mwSize n); + + + +/************************************************************************** + * Triangular back substitution. + * + * Solve the set of linear equations + * + * T*x = b + * + * where T is lower or upper triangular. + * + * Parameters: + * ul - 'U' for upper triangular, 'L' for lower triangular + * A - matrix of size n x m containing T + * b - vector of length k + * x - output vector of length k + * n - size of first dimension of A + * k - the size of the equation set, k<=n,m + * + * Note: + * The matrix A can be of any size n X m, as long as n,m >= k. + * Only the lower/upper triangle of the submatrix A(1:k,1:k) defines the + * matrix T (depending on the parameter ul). + * + **************************************************************************/ +void backsubst(char ul, double A[], double b[], double x[], mwSize n, mwSize k); + + + +/************************************************************************** + * Solve a set of equations using a Cholesky decomposition. + * + * Solve the set of linear equations + * + * M*x = b + * + * where M is positive definite with a known Cholesky decomposition: + * either M=L*L' (L lower triangular) or M=U'*U (U upper triangular). + * + * Parameters: + * ul - 'U' for upper triangular, 'L' for lower triangular decomposition + * A - matrix of size n x m with the Cholesky decomposition of M + * b - vector of length k + * x - output vector of length k + * n - size of first dimension of A + * k - the size of the equation set, k<=n,m + * + * Note: + * The matrix A can be of any size n X m, as long as n,m >= k. + * Only the lower/upper triangle of the submatrix A(1:k,1:k) is used as + * the Cholesky decomposition of M (depending on the parameter ul). + * + **************************************************************************/ +void cholsolve(char ul, double A[], double b[], double x[], mwSize n, mwSize k); + + + +/************************************************************************** + * Maximum absolute value. + * + * Returns the index of the coefficient with maximal absolute value in a vector. + * + * Parameters: + * x - vector of length n + * n - length of x + * + **************************************************************************/ +mwIndex maxabs(double x[], mwSize n); + + + +/************************************************************************** + * Maximum vector element. + * + * Returns the index of the maximal coefficient in a vector. + * + * Parameters: + * x - vector of length n + * n - length of x + * + **************************************************************************/ +mwIndex maxpos(double x[], mwSize n); + + + +/************************************************************************** + * Vector-vector dot product. + * + * Computes an operation of the form: + * + * c = a'*b + * + * Parameters: + * a, b - vectors of length n + * n - length of a,b + * + * Returns: The dot product c. + * + **************************************************************************/ +double dotprod(double a[], double b[], mwSize n); + + + +/************************************************************************** + * Indexed vector assignment. + * + * Perform a permutation assignment of the form + * + * y = x(ind) + * + * where ind is an array of indices to x. + * + * Parameters: + * y - output vector of length k + * x - input vector of arbitrary length + * ind - array of indices into x (indices begin at 0) + * k - length of the array ind + * + **************************************************************************/ +void vec_assign(double y[], double x[], mwIndex ind[], mwSize k); + + + +/************************************************************************** + * Matrix transpose. + * + * Computes Y := X' + * + * Parameters: + * X - input matrix of size n X m + * Y - output matrix of size m X n + * n, m - dimensions of X + * + **************************************************************************/ +void transpose(double X[], double Y[], mwSize n, mwSize m); + + + +/************************************************************************** + * Print a matrix. + * + * Parameters: + * A - matrix of size n X m + * n, m - dimensions of A + * matname - name of matrix to display + * + **************************************************************************/ +void printmat(double A[], int n, int m, char* matname); + + + +/************************************************************************** + * Print a sparse matrix. + * + * Parameters: + * A - sparse matrix of type double + * matname - name of matrix to display + * + **************************************************************************/ +void printspmat(mxArray *A, char* matname); + + +#endif +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/normcols.m Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,17 @@ +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)));
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/omp2mex.c Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,156 @@ +/************************************************************************** + * + * File name: omp2mex.c + * + * Ron Rubinstein + * Computer Science Department + * Technion, Haifa 32000 Israel + * ronrubin@cs + * + * Last Updated: 18.8.2009 + * + *************************************************************************/ + +#include "ompcore.h" +#include "omputils.h" +#include "mexutils.h" + + +/* Input Arguments */ + +#define IN_D prhs[0] +#define IN_X prhs[1] +#define IN_DtX prhs[2] +#define IN_XtX prhs[3] +#define IN_G prhs[4] +#define IN_EPS prhs[5] +#define IN_SPARSE_G prhs[6] +#define IN_MSGDELTA prhs[7] +#define IN_MAXATOMS prhs[8] +#define IN_PROFILE prhs[9] + + +/* Output Arguments */ + +#define GAMMA_OUT plhs[0] + + +/***************************************************************************************/ + + +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray*prhs[]) + +{ + double *D, *x, *DtX, *XtX, *G, eps, msgdelta; + int gmode, maxatoms, profile; + mwSize m, n, L; /* D is n x m , X is n x L, DtX is m x L */ + + + /* check parameters */ + + checkmatrix(IN_D, "OMP2", "D"); + checkmatrix(IN_X, "OMP2", "X"); + checkmatrix(IN_DtX, "OMP2", "DtX"); + checkmatrix(IN_XtX, "OMP2", "XtX"); + checkmatrix(IN_G, "OMP2", "G"); + + checkscalar(IN_EPS, "OMP2", "EPSILON"); + checkscalar(IN_SPARSE_G, "OMP2", "sparse_g"); + checkscalar(IN_MSGDELTA, "OMP2", "msgdelta"); + checkscalar(IN_MAXATOMS, "OMP2", "maxatoms"); + checkscalar(IN_PROFILE, "OMP2", "profile"); + + + /* get parameters */ + + x = D = DtX = XtX = G = 0; + + if (!mxIsEmpty(IN_D)) + D = mxGetPr(IN_D); + + if (!mxIsEmpty(IN_X)) + x = mxGetPr(IN_X); + + if (!mxIsEmpty(IN_DtX)) + DtX = mxGetPr(IN_DtX); + + if (!mxIsEmpty(IN_XtX)) + XtX = mxGetPr(IN_XtX); + + if (!mxIsEmpty(IN_G)) + G = mxGetPr(IN_G); + + eps = mxGetScalar(IN_EPS); + if ((int)(mxGetScalar(IN_SPARSE_G)+1e-2)) { + gmode = SPARSE_GAMMA; + } + else { + gmode = FULL_GAMMA; + } + msgdelta = mxGetScalar(IN_MSGDELTA); + if (mxGetScalar(IN_MAXATOMS) < -1e-5) { + maxatoms = -1; + } + else { + maxatoms = (int)(mxGetScalar(IN_MAXATOMS)+1e-2); + } + profile = (int)(mxGetScalar(IN_PROFILE)+1e-2); + + + /* check sizes */ + + if (D && x) { + n = mxGetM(IN_D); + m = mxGetN(IN_D); + L = mxGetN(IN_X); + + if (mxGetM(IN_X) != n) { + mexErrMsgTxt("D and X have incompatible sizes."); + } + + if (G) { + if (mxGetN(IN_G)!=mxGetM(IN_G)) { + mexErrMsgTxt("G must be a square matrix."); + } + if (mxGetN(IN_G) != m) { + mexErrMsgTxt("D and G have incompatible sizes."); + } + } + } + + else if (DtX && XtX) { + m = mxGetM(IN_DtX); + L = mxGetN(IN_DtX); + + /* set n to an arbitrary value that is at least the max possible number of selected atoms */ + + if (maxatoms>0) { + n = maxatoms; + } + else { + n = m; + } + + if ( !(mxGetM(IN_XtX)==L && mxGetN(IN_XtX)==1) && !(mxGetM(IN_XtX)==1 && mxGetN(IN_XtX)==L) ) { + mexErrMsgTxt("DtX and XtX have incompatible sizes."); + } + + if (mxGetN(IN_G)!=mxGetM(IN_G)) { + mexErrMsgTxt("G must be a square matrix."); + } + if (mxGetN(IN_G) != m) { + mexErrMsgTxt("DtX and G have incompatible sizes."); + } + } + + else { + mexErrMsgTxt("Either D and X, or DtX and XtX, must be specified."); + } + + + /* Do OMP! */ + + GAMMA_OUT = ompcore(D, x, DtX, XtX, G, n, m, L, maxatoms, eps, gmode, profile, msgdelta, 1); + + return; +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/omp2mex.m Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,23 @@ +%This is the Matlab interface to the OMP2 MEX implementation. +%The function is not for independent use, only through omp2.m. + + +%OMP2MEX Matlab interface to the OMP2 MEX implementation. +% GAMMA = OMP2MEX(D,X,DtX,XtX,G,EPSILON,SPARSE_G,MSGDELTA,MAXATOMS,PROFILE) +% invokes the OMP2 MEX function according to the specified parameters. Not +% all the parameters are required. Those among D, X, DtX, XtX and G which +% are not specified should be passed as []. +% +% EPSILON - the target error. +% SPARSE_G - returns a sparse GAMMA when nonzero, full GAMMA when zero. +% MSGDELTA - the delay in secs between messages. Zero means no messages. +% MAXATOMS - the max number of atoms per signal, negative for no max. +% PROFILE - nonzero means that profiling information should be printed. + + +% Ron Rubinstein +% Computer Science Department +% Technion, Haifa 32000 Israel +% ronrubin@cs +% +% April 2009
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/ompcore.c Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,409 @@ +/************************************************************************** + * + * File name: ompcore.c + * + * Ron Rubinstein + * Computer Science Department + * Technion, Haifa 32000 Israel + * ronrubin@cs + * + * Last Updated: 25.8.2009 + * + *************************************************************************/ + + +#include "ompcore.h" +#include "omputils.h" +#include "ompprof.h" +#include "myblas.h" +#include <math.h> +#include <string.h> + + + +/****************************************************************************** + * * + * Batch-OMP Implementation * + * * + ******************************************************************************/ + +mxArray* ompcore(double D[], double x[], double DtX[], double XtX[], double G[], mwSize n, mwSize m, mwSize L, + int T, double eps, int gamma_mode, int profile, double msg_delta, int erroromp) +{ + + profdata pd; + mxArray *Gamma; + mwIndex i, j, signum, pos, *ind, *gammaIr, *gammaJc, gamma_count; + mwSize allocated_coefs, allocated_cols; + int DtX_specified, XtX_specified, batchomp, standardomp, *selected_atoms; + double *alpha, *r, *Lchol, *c, *Gsub, *Dsub, sum, *gammaPr, *tempvec1, *tempvec2; + double eps2, resnorm, delta, deltaprev, secs_remain; + int mins_remain, hrs_remain; + clock_t lastprint_time, starttime; + + + + /*** status flags ***/ + + DtX_specified = (DtX!=0); /* indicates whether D'*x was provided */ + XtX_specified = (XtX!=0); /* indicates whether sum(x.*x) was provided */ + + standardomp = (G==0); /* batch-omp or standard omp are selected depending on availability of G */ + batchomp = !standardomp; + + + + /*** allocate output matrix ***/ + + + if (gamma_mode == FULL_GAMMA) { + + /* allocate full matrix of size m X L */ + + Gamma = mxCreateDoubleMatrix(m, L, mxREAL); + gammaPr = mxGetPr(Gamma); + gammaIr = 0; + gammaJc = 0; + } + else { + + /* allocate sparse matrix with room for allocated_coefs nonzeros */ + + /* for error-omp, begin with L*sqrt(n)/2 allocated nonzeros, otherwise allocate L*T nonzeros */ + allocated_coefs = erroromp ? (mwSize)(ceil(L*sqrt((double)n)/2.0) + 1.01) : L*T; + Gamma = mxCreateSparse(m, L, allocated_coefs, mxREAL); + gammaPr = mxGetPr(Gamma); + gammaIr = mxGetIr(Gamma); + gammaJc = mxGetJc(Gamma); + gamma_count = 0; + gammaJc[0] = 0; + } + + + /*** helper arrays ***/ + + alpha = (double*)mxMalloc(m*sizeof(double)); /* contains D'*residual */ + ind = (mwIndex*)mxMalloc(n*sizeof(mwIndex)); /* indices of selected atoms */ + selected_atoms = (int*)mxMalloc(m*sizeof(int)); /* binary array with 1's for selected atoms */ + c = (double*)mxMalloc(n*sizeof(double)); /* orthogonal projection result */ + + /* current number of columns in Dsub / Gsub / Lchol */ + allocated_cols = erroromp ? (mwSize)(ceil(sqrt((double)n)/2.0) + 1.01) : T; + + /* Cholesky decomposition of D_I'*D_I */ + Lchol = (double*)mxMalloc(n*allocated_cols*sizeof(double)); + + /* temporary vectors for various computations */ + tempvec1 = (double*)mxMalloc(m*sizeof(double)); + tempvec2 = (double*)mxMalloc(m*sizeof(double)); + + if (batchomp) { + /* matrix containing G(:,ind) - the columns of G corresponding to the selected atoms, in order of selection */ + Gsub = (double*)mxMalloc(m*allocated_cols*sizeof(double)); + } + else { + /* matrix containing D(:,ind) - the selected atoms from D, in order of selection */ + Dsub = (double*)mxMalloc(n*allocated_cols*sizeof(double)); + + /* stores the residual */ + r = (double*)mxMalloc(n*sizeof(double)); + } + + if (!DtX_specified) { + /* contains D'*x for the current signal */ + DtX = (double*)mxMalloc(m*sizeof(double)); + } + + + + /*** initializations for error omp ***/ + + if (erroromp) { + eps2 = eps*eps; /* compute eps^2 */ + if (T<0 || T>n) { /* unspecified max atom num - set max atoms to n */ + T = n; + } + } + + + + /*** initialize timers ***/ + + initprofdata(&pd); /* initialize profiling counters */ + starttime = clock(); /* record starting time for eta computations */ + lastprint_time = starttime; /* time of last status display */ + + + + /********************** perform omp for each signal **********************/ + + + + for (signum=0; signum<L; ++signum) { + + + /* initialize residual norm and deltaprev for error-omp */ + + if (erroromp) { + if (XtX_specified) { + resnorm = XtX[signum]; + } + else { + resnorm = dotprod(x+n*signum, x+n*signum, n); + addproftime(&pd, XtX_TIME); + } + deltaprev = 0; /* delta tracks the value of gamma'*G*gamma */ + } + else { + /* ignore residual norm stopping criterion */ + eps2 = 0; + resnorm = 1; + } + + + if (resnorm>eps2 && T>0) { + + /* compute DtX */ + + if (!DtX_specified) { + matT_vec(1, D, x+n*signum, DtX, n, m); + addproftime(&pd, DtX_TIME); + } + + + /* initialize alpha := DtX */ + + memcpy(alpha, DtX + m*signum*DtX_specified, m*sizeof(double)); + + + /* mark all atoms as unselected */ + + for (i=0; i<m; ++i) { + selected_atoms[i] = 0; + } + + } + + + /* main loop */ + + i=0; + while (resnorm>eps2 && i<T) { + + /* index of next atom */ + + pos = maxabs(alpha, m); + addproftime(&pd, MAXABS_TIME); + + + /* stop criterion: selected same atom twice, or inner product too small */ + + if (selected_atoms[pos] || alpha[pos]*alpha[pos]<1e-14) { + break; + } + + + /* mark selected atom */ + + ind[i] = pos; + selected_atoms[pos] = 1; + + + /* matrix reallocation */ + + if (erroromp && i>=allocated_cols) { + + allocated_cols = (mwSize)(ceil(allocated_cols*MAT_INC_FACTOR) + 1.01); + + Lchol = (double*)mxRealloc(Lchol,n*allocated_cols*sizeof(double)); + + batchomp ? (Gsub = (double*)mxRealloc(Gsub,m*allocated_cols*sizeof(double))) : + (Dsub = (double*)mxRealloc(Dsub,n*allocated_cols*sizeof(double))) ; + } + + + /* append column to Gsub or Dsub */ + + if (batchomp) { + memcpy(Gsub+i*m, G+pos*m, m*sizeof(double)); + } + else { + memcpy(Dsub+i*n, D+pos*n, n*sizeof(double)); + } + + + /*** Cholesky update ***/ + + if (i==0) { + *Lchol = 1; + } + else { + + /* incremental Cholesky decomposition: compute next row of Lchol */ + + if (standardomp) { + matT_vec(1, Dsub, D+n*pos, tempvec1, n, i); /* compute tempvec1 := Dsub'*d where d is new atom */ + addproftime(&pd, DtD_TIME); + } + else { + vec_assign(tempvec1, Gsub+i*m, ind, i); /* extract tempvec1 := Gsub(ind,i) */ + } + backsubst('L', Lchol, tempvec1, tempvec2, n, i); /* compute tempvec2 = Lchol \ tempvec1 */ + for (j=0; j<i; ++j) { /* write tempvec2 to end of Lchol */ + Lchol[j*n+i] = tempvec2[j]; + } + + /* compute Lchol(i,i) */ + sum = 0; + for (j=0; j<i; ++j) { /* compute sum of squares of last row without Lchol(i,i) */ + sum += SQR(Lchol[j*n+i]); + } + if ( (1-sum) <= 1e-14 ) { /* Lchol(i,i) is zero => selected atoms are dependent */ + break; + } + Lchol[i*n+i] = sqrt(1-sum); + } + + addproftime(&pd, LCHOL_TIME); + + i++; + + + /* perform orthogonal projection and compute sparse coefficients */ + + vec_assign(tempvec1, DtX + m*signum*DtX_specified, ind, i); /* extract tempvec1 = DtX(ind) */ + cholsolve('L', Lchol, tempvec1, c, n, i); /* solve LL'c = tempvec1 for c */ + addproftime(&pd, COMPCOEF_TIME); + + + /* update alpha = D'*residual */ + + if (standardomp) { + mat_vec(-1, Dsub, c, r, n, i); /* compute r := -Dsub*c */ + vec_sum(1, x+n*signum, r, n); /* compute r := x+r */ + + + /*memcpy(r, x+n*signum, n*sizeof(double)); /* assign r := x */ + /*mat_vec1(-1, Dsub, c, 1, r, n, i); /* compute r := r-Dsub*c */ + + addproftime(&pd, COMPRES_TIME); + matT_vec(1, D, r, alpha, n, m); /* compute alpha := D'*r */ + addproftime(&pd, DtR_TIME); + + /* update residual norm */ + if (erroromp) { + resnorm = dotprod(r, r, n); + addproftime(&pd, UPDATE_RESNORM_TIME); + } + } + else { + mat_vec(1, Gsub, c, tempvec1, m, i); /* compute tempvec1 := Gsub*c */ + memcpy(alpha, DtX + m*signum*DtX_specified, m*sizeof(double)); /* set alpha = D'*x */ + vec_sum(-1, tempvec1, alpha, m); /* compute alpha := alpha - tempvec1 */ + addproftime(&pd, UPDATE_DtR_TIME); + + /* update residual norm */ + if (erroromp) { + vec_assign(tempvec2, tempvec1, ind, i); /* assign tempvec2 := tempvec1(ind) */ + delta = dotprod(c,tempvec2,i); /* compute c'*tempvec2 */ + resnorm = resnorm - delta + deltaprev; /* residual norm update */ + deltaprev = delta; + addproftime(&pd, UPDATE_RESNORM_TIME); + } + } + } + + + /*** generate output vector gamma ***/ + + if (gamma_mode == FULL_GAMMA) { /* write the coefs in c to their correct positions in gamma */ + for (j=0; j<i; ++j) { + gammaPr[m*signum + ind[j]] = c[j]; + } + } + else { + /* sort the coefs by index before writing them to gamma */ + quicksort(ind,c,i); + addproftime(&pd, INDEXSORT_TIME); + + /* gamma is full - reallocate */ + if (gamma_count+i >= allocated_coefs) { + + while(gamma_count+i >= allocated_coefs) { + allocated_coefs = (mwSize)(ceil(GAMMA_INC_FACTOR*allocated_coefs) + 1.01); + } + + mxSetNzmax(Gamma, allocated_coefs); + mxSetPr(Gamma, mxRealloc(gammaPr, allocated_coefs*sizeof(double))); + mxSetIr(Gamma, mxRealloc(gammaIr, allocated_coefs*sizeof(mwIndex))); + + gammaPr = mxGetPr(Gamma); + gammaIr = mxGetIr(Gamma); + } + + /* append coefs to gamma and update the indices */ + for (j=0; j<i; ++j) { + gammaPr[gamma_count] = c[j]; + gammaIr[gamma_count] = ind[j]; + gamma_count++; + } + gammaJc[signum+1] = gammaJc[signum] + i; + } + + + + /*** display status messages ***/ + + if (msg_delta>0 && (clock()-lastprint_time)/(double)CLOCKS_PER_SEC >= msg_delta) + { + lastprint_time = clock(); + + /* estimated remainig time */ + secs2hms( ((L-signum-1)/(double)(signum+1)) * ((lastprint_time-starttime)/(double)CLOCKS_PER_SEC) , + &hrs_remain, &mins_remain, &secs_remain); + + mexPrintf("omp: signal %d / %d, estimated remaining time: %02d:%02d:%05.2f\n", + signum+1, L, hrs_remain, mins_remain, secs_remain); + mexEvalString("drawnow;"); + } + + } + + /* end omp */ + + + + /*** print final messages ***/ + + if (msg_delta>0) { + mexPrintf("omp: signal %d / %d\n", signum, L); + } + + if (profile) { + printprofinfo(&pd, erroromp, batchomp, L); + } + + + + /* free memory */ + + if (!DtX_specified) { + mxFree(DtX); + } + if (standardomp) { + mxFree(r); + mxFree(Dsub); + } + else { + mxFree(Gsub); + } + mxFree(tempvec2); + mxFree(tempvec1); + mxFree(Lchol); + mxFree(c); + mxFree(selected_atoms); + mxFree(ind); + mxFree(alpha); + + return Gamma; +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/ompcore.h Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,80 @@ +/************************************************************************** + * + * File name: ompcore.h + * + * Ron Rubinstein + * Computer Science Department + * Technion, Haifa 32000 Israel + * ronrubin@cs + * + * Last Updated: 18.8.2009 + * + * Contains the core implementation of Batch-OMP / OMP-Cholesky. + * + *************************************************************************/ + + +#ifndef __OMP_CORE_H__ +#define __OMP_CORE_H__ + + +#include "mex.h" + + + +/************************************************************************** + * Perform Batch-OMP or OMP-Cholesky on a specified set of signals, using + * either a fixed number of atoms or an error bound. + * + * Parameters (not all required): + * + * D - the dictionary, of size n X m + * x - the signals, of size n X L + * DtX - D'*x, of size m X L + * XtX - squared norms of the signals in x, sum(x.*x), of length L + * G - D'*D, of size m X m + * T - target sparsity, or maximal number of atoms for error-based OMP + * eps - target residual norm for error-based OMP + * gamma_mode - one of the constants FULL_GAMMA or SPARSE_GAMMA + * profile - if non-zero, profiling info is printed + * msg_delta - positive: the # of seconds between status prints, otherwise: nothing is printed + * erroromp - if nonzero indicates error-based OMP, otherwise fixed sparsity OMP + * + * Usage: + * + * The function can be called using different parameters, and will have + * different complexity depending on the parameters specified. Arrays which + * are not specified should be passed as null (0). When G is specified, + * Batch-OMP is performed. Otherwise, OMP-Cholesky is performed. + * + * Fixed-sparsity usage: + * --------------------- + * Either DtX, or D and x, must be specified. Specifying DtX is more efficient. + * XtX does not need to be specified. + * When D and x are specified, G is not required. However, not providing G + * will significantly degrade efficiency. + * The number of atoms must be specified in T. The value of eps is ignored. + * Finally, set erroromp to 0. + * + * Error-OMP usage: + * ---------------- + * Either DtX and Xtx, or D and x, must be specified. Specifying DtX and XtX + * is more efficient. + * When D and x are specified, G is not required. However, not providing G + * will significantly degrade efficiency. + * The target error must be specified in eps. A hard limit on the number + * of atoms can also be specified via the parameter T. Otherwise, T should + * be negative. Finally, set erroromp to nonzero. + * + * + * Returns: + * An mxArray containing the sparse representations of the signals in x + * (allocated using the appropriate mxCreateXXX() function). + * The array is either full or sparse, depending on gamma_mode. + * + **************************************************************************/ +mxArray* ompcore(double D[], double x[], double DtX[], double XtX[], double G[], mwSize n, mwSize m, mwSize L, + int T, double eps, int gamma_mode, int profile, double msg_delta, int erroromp); + + +#endif
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/ompmex.c Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,133 @@ +/************************************************************************** + * + * File name: ompmex.c + * + * Ron Rubinstein + * Computer Science Department + * Technion, Haifa 32000 Israel + * ronrubin@cs + * + * Last Updated: 18.8.2009 + * + *************************************************************************/ + +#include "ompcore.h" +#include "omputils.h" +#include "mexutils.h" + + +/* Input Arguments */ + +#define IN_D prhs[0] +#define IN_X prhs[1] +#define IN_DtX prhs[2] +#define IN_G prhs[3] +#define IN_T prhs[4] +#define IN_SPARSE_G prhs[5] +#define IN_MSGDELTA prhs[6] +#define IN_PROFILE prhs[7] + + +/* Output Arguments */ + +#define GAMMA_OUT plhs[0] + + +/***************************************************************************************/ + + +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray*prhs[]) + +{ + double *D, *x, *DtX, *G, msgdelta; + int gmode, profile, T; + mwSize m, n, L; /* D is n x m , X is n x L, DtX is m x L */ + + + /* check parameters */ + + checkmatrix(IN_D, "OMP", "D"); + checkmatrix(IN_X, "OMP", "X"); + checkmatrix(IN_DtX, "OMP", "DtX"); + checkmatrix(IN_G, "OMP", "G"); + + checkscalar(IN_T, "OMP", "T"); + checkscalar(IN_SPARSE_G, "OMP", "sparse_g"); + checkscalar(IN_MSGDELTA, "OMP", "msgdelta"); + checkscalar(IN_PROFILE, "OMP", "profile"); + + + /* get parameters */ + + x = D = DtX = G = 0; + + if (!mxIsEmpty(IN_D)) + D = mxGetPr(IN_D); + + if (!mxIsEmpty(IN_X)) + x = mxGetPr(IN_X); + + if (!mxIsEmpty(IN_DtX)) + DtX = mxGetPr(IN_DtX); + + if (!mxIsEmpty(IN_G)) + G = mxGetPr(IN_G); + + T = (int)(mxGetScalar(IN_T)+1e-2); + if ((int)(mxGetScalar(IN_SPARSE_G)+1e-2)) { + gmode = SPARSE_GAMMA; + } + else { + gmode = FULL_GAMMA; + } + msgdelta = mxGetScalar(IN_MSGDELTA); + profile = (int)(mxGetScalar(IN_PROFILE)+1e-2); + + + /* check sizes */ + + if (D && x) { + n = mxGetM(IN_D); + m = mxGetN(IN_D); + L = mxGetN(IN_X); + + if (mxGetM(IN_X) != n) { + mexErrMsgTxt("D and X have incompatible sizes."); + } + + if (G) { + if (mxGetN(IN_G)!=mxGetM(IN_G)) { + mexErrMsgTxt("G must be a square matrix."); + } + if (mxGetN(IN_G) != m) { + mexErrMsgTxt("D and G have incompatible sizes."); + } + } + } + + else if (DtX) { + m = mxGetM(IN_DtX); + L = mxGetN(IN_DtX); + + n = T; /* arbitrary - it is enough to assume signal length is T */ + + if (mxGetN(IN_G)!=mxGetM(IN_G)) { + mexErrMsgTxt("G must be a square matrix."); + } + if (mxGetN(IN_G) != m) { + mexErrMsgTxt("DtX and G have incompatible sizes."); + } + } + + else { + mexErrMsgTxt("Either D and X, or DtX, must be specified."); + } + + + /* Do OMP! */ + + GAMMA_OUT = ompcore(D, x, DtX, 0, G, n, m, L, T, 0, gmode, profile, msgdelta, 0); + + return; +} +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/ompmex.m Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,22 @@ +%This is the Matlab interface to the OMP MEX implementation. +%The function is not for independent use, only through omp.m. + + +%OMPMEX Matlab interface to the OMP MEX implementation. +% GAMMA = OMPMEX(D,X,DtX,G,L,SPARSE_G,MSGDELTA,PROFILE) invokes the OMP +% MEX function according to the specified parameters. Not all the +% parameters are required. Those among D, X, DtX and G which are not +% specified should be passed as []. +% +% L - the target sparsity. +% SPARSE_G - returns a sparse GAMMA when nonzero, full GAMMA when zero. +% MSGDELTA - the delay in secs between messages. Zero means no messages. +% PROFILE - nonzero means that profiling information should be printed. + + +% Ron Rubinstein +% Computer Science Department +% Technion, Haifa 32000 Israel +% ronrubin@cs +% +% April 2009
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/ompprof.c Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,113 @@ +/************************************************************************** + * + * File name: ompprof.c + * + * Ron Rubinstein + * Computer Science Department + * Technion, Haifa 32000 Israel + * ronrubin@cs + * + * Last Updated: 11.4.2009 + * + *************************************************************************/ + + +#include "ompprof.h" + + +/* initialize profiling information */ + +void initprofdata(profdata *pd) +{ + pd->DtX_time = 0; + pd->XtX_time = 0; + pd->DtR_time = 0; + pd->maxabs_time = 0; + pd->DtD_time = 0; + pd->Lchol_time = 0; + pd->compcoef_time = 0; + pd->update_DtR_time = 0; + pd->update_resnorm_time = 0; + pd->compres_time = 0; + pd->indexsort_time = 0; + + pd->DtX_time_counted = 0; + pd->XtX_time_counted = 0; + pd->DtR_time_counted = 0; + pd->DtD_time_counted = 0; + pd->update_DtR_time_counted = 0; + pd->resnorm_time_counted = 0; + pd->compres_time_counted = 0; + pd->indexsort_time_counted = 0; + + pd->prevtime = clock(); +} + + +/* add elapsed time to profiling data according to specified computation */ + +void addproftime(profdata *pd, int comptype) +{ + switch(comptype) { + case DtX_TIME: pd->DtX_time += clock()-pd->prevtime; pd->DtX_time_counted = 1; break; + case XtX_TIME: pd->XtX_time += clock()-pd->prevtime; pd->XtX_time_counted = 1; break; + case DtR_TIME: pd->DtR_time += clock()-pd->prevtime; pd->DtR_time_counted = 1; break; + case DtD_TIME: pd->DtD_time += clock()-pd->prevtime; pd->DtD_time_counted = 1; break; + case COMPRES_TIME: pd->compres_time += clock()-pd->prevtime; pd->compres_time_counted = 1; break; + case UPDATE_DtR_TIME: pd->update_DtR_time += clock()-pd->prevtime; pd->update_DtR_time_counted = 1; break; + case UPDATE_RESNORM_TIME: pd->update_resnorm_time += clock()-pd->prevtime; pd->resnorm_time_counted = 1; break; + case INDEXSORT_TIME: pd->indexsort_time += clock()-pd->prevtime; pd->indexsort_time_counted = 1; break; + case MAXABS_TIME: pd->maxabs_time += clock()-pd->prevtime; break; + case LCHOL_TIME: pd->Lchol_time += clock()-pd->prevtime; break; + case COMPCOEF_TIME: pd->compcoef_time += clock()-pd->prevtime; break; + } + pd->prevtime = clock(); +} + + +/* print profiling info */ + +void printprofinfo(profdata *pd, int erroromp, int batchomp, int signum) +{ + clock_t tottime; + + tottime = pd->DtX_time + pd->XtX_time + pd->DtR_time + pd->DtD_time + pd->compres_time + pd->maxabs_time + + pd->Lchol_time + pd->compcoef_time + pd->update_DtR_time + pd->update_resnorm_time + pd->indexsort_time; + + mexPrintf("\n\n***** Profiling information for %s *****\n\n", erroromp? "OMP2" : "OMP"); + + mexPrintf("OMP mode: %s\n\n", batchomp? "Batch-OMP" : "OMP-Cholesky"); + + mexPrintf("Total signals processed: %d\n\n", signum); + + if (pd->DtX_time_counted) { + mexPrintf("Compute DtX time: %7.3lf seconds\n", pd->DtX_time/(double)CLOCKS_PER_SEC); + } + if (pd->XtX_time_counted) { + mexPrintf("Compute XtX time: %7.3lf seconds\n", pd->XtX_time/(double)CLOCKS_PER_SEC); + } + mexPrintf("Max abs time: %7.3lf seconds\n", pd->maxabs_time/(double)CLOCKS_PER_SEC); + if (pd->DtD_time_counted) { + mexPrintf("Compute DtD time: %7.3lf seconds\n", pd->DtD_time/(double)CLOCKS_PER_SEC); + } + mexPrintf("Lchol update time: %7.3lf seconds\n", pd->Lchol_time/(double)CLOCKS_PER_SEC); + mexPrintf("Compute coef time: %7.3lf seconds\n", pd->compcoef_time/(double)CLOCKS_PER_SEC); + if (pd->compres_time_counted) { + mexPrintf("Compute R time: %7.3lf seconds\n", pd->compres_time/(double)CLOCKS_PER_SEC); + } + if (pd->DtR_time_counted) { + mexPrintf("Compute DtR time: %7.3lf seconds\n", pd->DtR_time/(double)CLOCKS_PER_SEC); + } + if (pd->update_DtR_time_counted) { + mexPrintf("Update DtR time: %7.3lf seconds\n", pd->update_DtR_time/(double)CLOCKS_PER_SEC); + } + if (pd->resnorm_time_counted) { + mexPrintf("Update resnorm time: %7.3lf seconds\n", pd->update_resnorm_time/(double)CLOCKS_PER_SEC); + } + if (pd->indexsort_time_counted) { + mexPrintf("Index sort time: %7.3lf seconds\n", pd->indexsort_time/(double)CLOCKS_PER_SEC); + } + mexPrintf("---------------------------------------\n"); + mexPrintf("Total time: %7.3lf seconds\n\n", tottime/(double)CLOCKS_PER_SEC); +} +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/ompprof.h Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,106 @@ +/************************************************************************** + * + * File name: ompprof.h + * + * Ron Rubinstein + * Computer Science Department + * Technion, Haifa 32000 Israel + * ronrubin@cs + * + * Last Updated: 18.8.2009 + * + * Collection of definitions and functions for profiling the OMP method. + * + *************************************************************************/ + + +#ifndef __OMP_PROF_H__ +#define __OMP_PROF_H__ + +#include "mex.h" +#include <time.h> + + + +/************************************************************************** + * + * Constants and data types. + * + **************************************************************************/ + + +/* constants denoting the various parts of the algorithm */ + +enum { DtX_TIME, XtX_TIME, DtR_TIME, MAXABS_TIME, DtD_TIME, LCHOL_TIME, COMPCOEF_TIME, + UPDATE_DtR_TIME, UPDATE_RESNORM_TIME, COMPRES_TIME, INDEXSORT_TIME }; + + + +/* profiling data container with counters for each part of the algorithm */ + +typedef struct profdata +{ + clock_t prevtime; /* the time when last initialization/call to addproftime() was performed */ + + clock_t DtX_time; + clock_t XtX_time; + clock_t DtR_time; + clock_t maxabs_time; + clock_t DtD_time; + clock_t Lchol_time; + clock_t compcoef_time; + clock_t update_DtR_time; + clock_t update_resnorm_time; + clock_t compres_time; + clock_t indexsort_time; + + /* flags indicating whether profiling data was gathered */ + int DtX_time_counted; + int XtX_time_counted; + int DtR_time_counted; + int DtD_time_counted; + int update_DtR_time_counted; + int resnorm_time_counted; + int compres_time_counted; + int indexsort_time_counted; + +} profdata; + + + +/************************************************************************** + * + * Initialize a profdata structure, zero all counters, and start its timer. + * + **************************************************************************/ +void initprofdata(profdata *pd); + + +/************************************************************************** + * + * Add elapsed time from last call to addproftime(), or from initialization + * of profdata, to the counter specified by comptype. comptype must be one + * of the constants in the enumeration above. + * + **************************************************************************/ +void addproftime(profdata *pd, int comptype); + + +/************************************************************************** + * + * Print the current contents of the counters in profdata. + * + * Parameters: + * pd - the profdata to print + * erroromp - indicates whether error-based (nonzero) or sparsity-based (zero) + * omp was performed. + * batchomp - indicates whether batch-omp (nonzero) or omp-cholesky (zero) + * omp was performed. + * signum - number of signals processed by omp + * + **************************************************************************/ +void printprofinfo(profdata *pd, int erroromp, int batchomp, int signum); + + +#endif +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/omputils.c Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,89 @@ +/************************************************************************** + * + * File name: omputils.c + * + * Ron Rubinstein + * Computer Science Department + * Technion, Haifa 32000 Israel + * ronrubin@cs + * + * Last Updated: 18.8.2009 + * + *************************************************************************/ + +#include "omputils.h" +#include <math.h> + + +const char FULL_GAMMA_STR[] = "full"; +const char SPARSE_GAMMA_STR[] = "sparse"; + + +/* convert seconds to hours, minutes and seconds */ + +void secs2hms(double sectot, int *hrs, int *mins, double *secs) +{ + *hrs = (int)(floor(sectot/3600)+1e-2); + sectot = sectot - 3600*(*hrs); + *mins = (int)(floor(sectot/60)+1e-2); + *secs = sectot - 60*(*mins); +} + + +/* quicksort, public-domain C implementation by Darel Rex Finley. */ +/* modification: sorts the array data[] as well, according to the values in the array vals[] */ + +#define MAX_LEVELS 300 + +void quicksort(mwIndex vals[], double data[], mwIndex n) { + + long piv, beg[MAX_LEVELS], end[MAX_LEVELS], i=0, L, R, swap ; + double datapiv; + + beg[0]=0; + end[0]=n; + + while (i>=0) { + + L=beg[i]; + R=end[i]-1; + + if (L<R) { + + piv=vals[L]; + datapiv=data[L]; + + while (L<R) + { + while (vals[R]>=piv && L<R) + R--; + if (L<R) { + vals[L]=vals[R]; + data[L++]=data[R]; + } + + while (vals[L]<=piv && L<R) + L++; + if (L<R) { + vals[R]=vals[L]; + data[R--]=data[L]; + } + } + + vals[L]=piv; + data[L]=datapiv; + + beg[i+1]=L+1; + end[i+1]=end[i]; + end[i++]=L; + + if (end[i]-beg[i] > end[i-1]-beg[i-1]) { + swap=beg[i]; beg[i]=beg[i-1]; beg[i-1]=swap; + swap=end[i]; end[i]=end[i-1]; end[i-1]=swap; + } + } + else { + i--; + } + } +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/omputils.h Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,77 @@ +/************************************************************************** + * + * File name: omputils.h + * + * Ron Rubinstein + * Computer Science Department + * Technion, Haifa 32000 Israel + * ronrubin@cs + * + * Last Updated: 18.8.2009 + * + * Utility definitions and functions for the OMP library. + * + *************************************************************************/ + + +#ifndef __OMP_UTILS_H__ +#define __OMP_UTILS_H__ + +#include "mex.h" + + +/* constants for the representation mode of gamma */ + +extern const char FULL_GAMMA_STR[]; /* "full" */ +extern const char SPARSE_GAMMA_STR[]; /* "sparse" */ + + +#define FULL_GAMMA 1 +#define SPARSE_GAMMA 2 +#define INVALID_MODE 3 + + + +/************************************************************************** + * Memory management for OMP2. + * + * GAMMA_INC_FACTOR: + * The matrix GAMMA is allocated with sqrt(n)/2 coefficients per signal, + * for a total of nzmax = L*sqrt(n)/2 nonzeros. Whenever GAMMA needs to be + * increased, it is increased by a factor of GAMMA_INC_FACTOR. + * + * MAT_INC_FACTOR: + * The matrices Lchol, Gsub and Dsub are allocated with sqrt(n)/2 + * columns each. If additional columns are needed, this number is + * increased by a factor of MAT_INC_FACTOR. + **************************************************************************/ + +#define GAMMA_INC_FACTOR (1.4) +#define MAT_INC_FACTOR (1.6) + + + +/************************************************************************** + * Convert number of seconds to hour, minute and second representation. + * + * Parameters: + * sectot - total number of seconds + * hrs, mins, secs - output hours (whole) and minutes (whole) and seconds + * + **************************************************************************/ +void secs2hms(double sectot, int *hrs, int *mins, double *secs); + + + +/************************************************************************** + * QuickSort - public-domain C implementation by Darel Rex Finley. + * + * Modified to sort both the array vals[] and the array data[] according + * to the values in the array vals[]. + * + **************************************************************************/ +void quicksort(mwIndex vals[], double data[], mwIndex n); + + +#endif +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/printf.m Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,26 @@ +function str = printf(varargin) +%PRINTF Print formatted text to screen. +% PRINTF(FMT,VAL1,VAL2,...) formats the data in VAL1,VAL2,... according to +% the format string FMT, and prints the result to the screen. +% +% The call to PRINTF(FMT,VAL1,VAL2,...) simply invokes the call +% DISP(SPRINTF(FMT,VAL1,VAL2,...)). For a complete description of the +% format string options see function SPRINTF. +% +% STR = PRINTF(...) also returns the formatted string. + + +% Ron Rubinstein +% Computer Science Department +% Technion, Haifa 32000 Israel +% ronrubin@cs +% +% April 2008 + + +if (nargout>0) + str = sprintf(varargin{:}); + disp(str); +else + disp(sprintf(varargin{:})); +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/private/make.m Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,41 @@ +function make +%MAKE Build the OMPBox package. +% MAKE compiles all OMPBox MEX functions, using Matlab's default MEX +% compiler. If the MEX compiler has not been set-up before, please run +% +% mex -setup +% +% before using this MAKE file. + +% Ron Rubinstein +% Computer Science Department +% Technion, Haifa 32000 Israel +% ronrubin@cs +% +% August 2009 + + +% detect platform + +compstr = computer; +is64bit = strcmp(compstr(end-1:end),'64'); + + +% compilation parameters + +compile_params = cell(0); +if (is64bit) + compile_params{1} = '-largeArrayDims'; +end + + +% Compile files % + +ompsources = {'mexutils.c','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{:}); +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/private/mexutils.c Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,79 @@ +/************************************************************************** + * + * File name: mexutils.c + * + * Ron Rubinstein + * Computer Science Department + * Technion, Haifa 32000 Israel + * ronrubin@cs + * + * Last Updated: 15.8.2009 + * + *************************************************************************/ + +#include "mexutils.h" +#include <math.h> + + + +/* verify that the mxArray contains a double matrix */ + +void checkmatrix(const mxArray *param, char *fname, char *pname) +{ + char errmsg[100]; + sprintf(errmsg, "%.15s requires that %.25s be a double matrix.", fname, pname); + if (!mxIsDouble(param) || mxIsComplex(param) || mxGetNumberOfDimensions(param)>2) { + mexErrMsgTxt(errmsg); + } +} + + +/* verify that the mxArray contains a 1-D double vector */ + +void checkvector(const mxArray *param, char *fname, char *pname) +{ + char errmsg[100]; + sprintf(errmsg, "%.15s requires that %.25s be a double vector.", fname, pname); + if (!mxIsDouble(param) || mxIsComplex(param) || mxGetNumberOfDimensions(param)>2 || (mxGetM(param)!=1 && mxGetN(param)!=1)) { + mexErrMsgTxt(errmsg); + } +} + + +/* verify that the mxArray contains a double scalar */ + +void checkscalar(const mxArray *param, char *fname, char *pname) +{ + char errmsg[100]; + sprintf(errmsg, "%.15s requires that %.25s be a double scalar.", fname, pname); + if (!mxIsDouble(param) || mxIsComplex(param) || mxGetNumberOfDimensions(param)>2 || + mxGetM(param)!=1 || mxGetN(param)!=1) + { + mexErrMsgTxt(errmsg); + } +} + + +/* verify that the mxArray contains a sparse matrix */ + +void checksparse(const mxArray *param, char *fname, char *pname) +{ + char errmsg[100]; + sprintf(errmsg, "%.15s requires that %.25s be sparse.", fname, pname); + if (!mxIsSparse(param)) { + mexErrMsgTxt(errmsg); + } +} + + +/* verify that the mxArray contains a 1-dimensional cell array */ + +void checkcell_1d(const mxArray *param, char *fname, char *pname) +{ + char errmsg[100]; + sprintf(errmsg, "%.15s requires that %.25s be a 1-D cell array.", fname, pname); + if (!mxIsCell(param) || mxGetNumberOfDimensions(param)>2 || (mxGetM(param)!=1 && mxGetN(param)!=1)) { + mexErrMsgTxt(errmsg); + } +} +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/private/mexutils.h Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,103 @@ +/************************************************************************** + * + * File name: mexutils.h + * + * Ron Rubinstein + * Computer Science Department + * Technion, Haifa 32000 Israel + * ronrubin@cs + * + * Last Updated: 18.8.2009 + * + * Utility functions for MEX files. + * + *************************************************************************/ + + +#ifndef __MEX_UTILS_H__ +#define __MEX_UTILS_H__ + +#include "mex.h" + + + +/************************************************************************** + * Function checkmatrix: + * + * Verify that the specified mxArray is real, of type double, and has + * no more than two dimensions. If not, an error message is printed + * and the mex file terminates. + * + * Parameters: + * param - the mxArray to be checked + * fname - the name of the function where the error occured (15 characters or less) + * pname - the name of the parameter (25 characters or less) + * + **************************************************************************/ +void checkmatrix(const mxArray *param, char *fname, char *pname); + + +/************************************************************************** + * Function checkvector: + * + * Verify that the specified mxArray is 1-D, real, and of type double. The + * vector may be a column or row vector. Otherwise, an error message is + * printed and the mex file terminates. + * + * Parameters: + * param - the mxArray to be checked + * fname - the name of the function where the error occured (15 characters or less) + * pname - the name of the parameter (25 characters or less) + * + **************************************************************************/ +void checkvector(const mxArray *param, char *fname, char *pname); + + +/************************************************************************** + * Function checkscalar: + * + * Verify that the specified mxArray represents a real double scalar value. + * If not, an error message is printed and the mex file terminates. + * + * Parameters: + * param - the mxArray to be checked + * fname - the name of the function where the error occured (15 characters or less) + * pname - the name of the parameter (25 characters or less) + * + **************************************************************************/ +void checkscalar(const mxArray *param, char *fname, char *pname); + + +/************************************************************************** + * Function checksparse: + * + * Verify that the specified mxArray contains a sparse matrix. If not, + * an error message is printed and the mex file terminates. + * + * Parameters: + * param - the mxArray to be checked + * fname - the name of the function where the error occured (15 characters or less) + * pname - the name of the parameter (25 characters or less) + * + **************************************************************************/ +void checksparse(const mxArray *param, char *fname, char *pname); + + +/************************************************************************** + * Function checkcell_1d: + * + * Verify that the specified mxArray is a 1-D cell array. The cell array + * may be arranged as either a column or a row. If not, an error message + * is printed and the mex file terminates. + * + * Parameters: + * param - the mxArray to be checked + * fname - the name of the function where the error occured (15 characters or less) + * pname - the name of the parameter (25 characters or less) + * + **************************************************************************/ +void checkcell_1d(const mxArray *param, char *fname, char *pname); + + +#endif +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/private/printf.m Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,26 @@ +function str = printf(varargin) +%PRINTF Print formatted text to screen. +% PRINTF(FMT,VAL1,VAL2,...) formats the data in VAL1,VAL2,... according to +% the format string FMT, and prints the result to the screen. +% +% The call to PRINTF(FMT,VAL1,VAL2,...) simply invokes the call +% DISP(SPRINTF(FMT,VAL1,VAL2,...)). For a complete description of the +% format string options see function SPRINTF. +% +% STR = PRINTF(...) also returns the formatted string. + + +% Ron Rubinstein +% Computer Science Department +% Technion, Haifa 32000 Israel +% ronrubin@cs +% +% April 2008 + + +if (nargout>0) + str = sprintf(varargin{:}); + disp(str); +else + disp(sprintf(varargin{:})); +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/pwsmoothfield.m Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,32 @@ +function f = pwsmoothfield(n,var,alpha) +% PWSMOOTHFIELD(N,VAR,ALPHA) +% Generate an image of piecewise smooth filtered white noise +% +% N = sidelength of the field in pixels +% VAR = variance of original white nose +% ALPHA = fraction of FFT coefficents to keep +% +% Returns an N-by-N array +% +% This file is used with the kind permission of Stephen J. Wright +% (swright@cs.wisc.edu), and was originally included in the GPSR +% v1.0 distribution: http://www.lx.it.pt/~mtf/GPSR . + +% $Id: pwsmoothfield.m 1040 2008-06-26 20:29:02Z ewout78 $ + +f = sqrt(var)*randn(n); +F = fft2(f); +a = ceil(n*alpha/2); +b = fix(n*(1-alpha)); +F(a+1:a+b,:) = 0; +F(:,a+1:a+b) = 0; +f = real(ifft2(F)); + +for i = 1:n + for j = 1:n + if (j/n >= 15*(i/n - 0.5)^3 + 0.4) + f(i,j) = f(i,j) + sqrt(var); + end + end +end +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/reggrid.m Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,136 @@ +function [varargout] = reggrid(sz,num,mode) +%REGGRID Regular sampling grid. +% [I1,I2,...,Ip] = REGGRID([N1 N2 ... Np], NUM) returns the indices +% of a regular uniform sampling grid over a p-dimensional matrix with +% dimensions N1xN2x...xNp. NUM is the minimal number of required samples, +% and it is ensured that the actual number of samples, given by +% length(I1)xlength(I2)x...xlength(Ip), is at least as large as NUM. +% +% [I1,I2,...,Ip] = REGGRID([N1 N2 ... Np], NUM,'MODE') specifies the +% method for distributing the samples along each dimension. Valid modes +% include 'eqdist' (the default mode) and 'eqnum'. 'eqdist' indicates an +% equal distance between the samples in each dimension, while 'eqnum' +% indicates an equal number of samples in each dimension. +% +% Notes about MODE: +% +% 1. The 'eqnum' mode will generally fail when the p-th root of NUM +% (i.e. NUM^(1/p)) is larger than min([N1 N2 ... Np]). Thus 'eqdist' is +% the more useful choice for sampling an arbitrary number of samples +% from the matrix (up to the total number of matrix entries). +% +% 2. In both modes, the equality (of the distance between samples, or +% the number of samples in each dimension) is only approximate. This is +% because REGGRID attempts to maintain the appropriate equality while at +% the same time find a sampling pattern where the total number of +% samples is as close as possible to NUM. In general, the larger {Ni} +% and NUM are, the tighter the equality. +% +% Example: Sample a set of blocks uniformly from a 2D image. +% +% n = 512; blocknum = 20000; blocksize = [8 8]; +% im = rand(n,n); +% [i1,i2] = reggrid(size(im)-blocksize+1, blocknum); +% blocks = sampgrid(im, blocksize, i1, i2); +% +% See also SAMPGRID. + +% Ron Rubinstein +% Computer Science Department +% Technion, Haifa 32000 Israel +% ronrubin@cs +% +% November 2007 + +dim = length(sz); + +if (nargin<3) + mode = 'eqdist'; +end + +if (any(sz<1)) + error(['Invalid matrix size : [' num2str(sz) ']']); +end + +if (num > prod(sz)) + warning(['Invalid number of samples, returning maximum number of samples.']); +elseif (num <= 0) + if (num < 0) + warning('Invalid number of samples, assuming 0 samples.'); + end + for i = 1:length(sz) + varargout{i} = []; + end + return; +end + + +if (strcmp(mode,'eqdist')) + + % approximate distance between samples: total volume divided by number of + % samples gives the average volume per sample. then, taking the p-th root + % gives the average distance between samples + d = (prod(sz)/num)^(1/dim); + + % compute the initial guess for number of samples in each dimension. + % then, while total number of samples is too large, decrese the number of + % samples by one in the dimension where the samples are the most crowded. + % finally, do the opposite process until just passing num, so the final + % number of samples is the closest to num from above. + + n = min(max(round(sz/d),1),sz); % set n so that it saturates at 1 and sz + + active_dims = find(n>1); % dimensions where the sample num can be reduced + while(prod(n)>num && ~isempty(active_dims)) + [y,id] = min((sz(active_dims)-1)./n(active_dims)); + n(active_dims(id)) = n(active_dims(id))-1; + if (n(active_dims(id)) < 2) + active_dims = find(n>1); + end + end + + active_dims = find(n<sz); % dimensions where the sample num can be increased + while(prod(n)<num && ~isempty(active_dims)) + [y,id] = max((sz(active_dims)-1)./n(active_dims)); + n(active_dims(id)) = n(active_dims(id))+1; + if (n(active_dims(id)) >= sz(active_dims(id))) + active_dims = find(n<sz); + end + end + + for i = 1:dim + varargout{i} = round((1:n(i))/n(i)*sz(i)); + varargout{i} = varargout{i} - floor((varargout{i}(1)-1)/2); + end + +elseif (strcmp(mode,'eqnum')) + + % same idea as above + n = min(max( ones(size(sz)) * round(num^(1/dim)) ,1),sz); + + active_dims = find(n>1); + while(prod(n)>num && ~isempty(active_dims)) + [y,id] = min((sz(active_dims)-1)./n(active_dims)); + n(active_dims(id)) = n(active_dims(id))-1; + if (n(active_dims(id)) < 2) + active_dims = find(n>1); + end + end + + active_dims = find(n<sz); + while(prod(n)<num && ~isempty(active_dims)) + [y,id] = max((sz(active_dims)-1)./n(active_dims)); + n(active_dims(id)) = n(active_dims(id))+1; + if (n(active_dims(id)) >= sz(active_dims(id))) + active_dims = find(n<sz); + end + end + + for i = 1:dim + varargout{i} = round((1:n(i))/n(i)*sz(i)); + varargout{i} = varargout{i} - floor((varargout{i}(1)-1)/2); + end +else + error('Invalid sampling mode'); +end +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/remove_dc.m Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,31 @@ +function [y,dc] = remove_dc(x,columns) +%REMOVE_DC Remove DC channel from signals. +% [Y,DC] = REMOVE_DC(X) removes the DC channel (i.e. the mean) from the +% specified (possibly multi-dimensional) signal X. Y is the DC-free +% signal and is the same size as X. DC is a scalar containing the mean of +% the signal X. +% +% [Y,DC] = REMOVE_DC(X,'columns') where X is a 2D matrix, treats the +% columns of X as a set of 1D signals, removing the DC channel from each +% one individually. Y is the same size as X and contains the DC-free +% signals. DC is a row vector of length size(X,2) containing the means of +% the signals in X. +% +% See also ADD_DC. + + +if (nargin==2 && strcmpi(columns,'columns')), columns = 1; +else columns = 0; +end + +if (columns) + dc = mean(x); + y = addtocols(x,-dc); +else + if (ndims(x)==2) % temporary, will remove in future + warning('Treating 2D matrix X as a single signal and not each column individually'); + end + dc = mean(x(:)); + y = x-dc; +end +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/rowlincomb.c Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,148 @@ +/************************************************************************** + * + * File name: rowlincomb.c + * + * Ron Rubinstein + * Computer Science Department + * Technion, Haifa 32000 Israel + * ronrubin@cs + * + * Last Updated: 21.5.2009 + * + *************************************************************************/ + +#include "mex.h" + + +/* Input Arguments */ + +#define X_IN prhs[0] +#define A_IN prhs[1] +#define ROWS_IN prhs[2] +#define COLS_IN prhs[3] + + +/* Output Arguments */ + +#define Y_OUT plhs[0] + + +void mexFunction(int nlhs, mxArray *plhs[], +int nrhs, const mxArray*prhs[]) + +{ + double *A, *x, *y, *rows, *cols; + mwSize m,n,m1,n1,m2,n2,rownum,colnum; + mwIndex *row_ids,*col_ids,i,j; + int colnumspecified=0; + + + /* Check for proper number of arguments */ + + if (nrhs!=3 && nrhs!=4) { + mexErrMsgTxt("Invalid number of input arguments."); + } else if (nlhs > 1) { + mexErrMsgTxt("Too many output arguments."); + } + + + /* Check the input dimensions */ + + m = mxGetM(A_IN); + n = mxGetN(A_IN); + if (!mxIsDouble(A_IN) || mxIsComplex(A_IN) || mxGetNumberOfDimensions(A_IN)>2) { + mexErrMsgTxt("ROWLINCOMB requires that A be a double matrix."); + } + + m1 = mxGetM(ROWS_IN); + n1 = mxGetN(ROWS_IN); + if (!mxIsDouble(ROWS_IN) || mxIsComplex(ROWS_IN) || (m1!=1 && n1!=1)) { + mexErrMsgTxt("ROWLINCOMB requires that ROWS be an index vector of type double."); + } + rownum = (m1 > n1) ? m1 : n1; /* the number of rows in the linear combination */ + + m2 = mxGetM(X_IN); + n2 = mxGetN(X_IN); + if (!mxIsDouble(X_IN) || mxIsComplex(X_IN) || ((m2!=1) && (n2!=1))) { + mexErrMsgTxt("ROWLINCOMB requires that X be a double vector."); + } + + if (m2 != rownum && n2 != rownum) { + mexErrMsgTxt("The length of X does not match the number of rows in ROWS."); + } + + if (nrhs==4) { + m1 = mxGetM(COLS_IN); + n1 = mxGetN(COLS_IN); + if (!mxIsDouble(COLS_IN) || mxIsComplex(COLS_IN) || (m1!=1 && n1!=1)) { + mexErrMsgTxt("ROWLINCOMB requires that COLS be an index vector of type double."); + } + colnum = (m1 > n1) ? m1 : n1; /* the number of columns */ + colnumspecified = 1; + cols = mxGetPr(COLS_IN); + + Y_OUT = mxCreateDoubleMatrix(1, colnum, mxREAL); + } + else { + cols = 0; + Y_OUT = mxCreateDoubleMatrix(1, n, mxREAL); + } + + + /* Assign pointers to the various parameters */ + A = mxGetPr(A_IN); + rows = mxGetPr(ROWS_IN); + x = mxGetPr(X_IN); + y = mxGetPr(Y_OUT); + + + /* check row indices */ + + row_ids = (mwIndex*)mxMalloc(rownum*sizeof(mwIndex)); + + for (i=0; i<rownum; ++i) { + row_ids[i] = (mwIndex)(rows[i]+0.1)-1; + if (row_ids[i]<0 || row_ids[i]>=m) { + mexErrMsgTxt("Row index in ROWS is out of range."); + } + } + + + + if (colnumspecified) { + + /* check column indices */ + col_ids = (mwIndex*)mxMalloc(colnum*sizeof(mwIndex)); + + for (i=0; i<colnum; ++i) { + col_ids[i] = (mwIndex)(cols[i]+0.1)-1; + if (col_ids[i]<0 || col_ids[i]>=n) { + mexErrMsgTxt("Column index in COLS is out of range."); + } + } + + /* Do the actual computation */ + for (j=0; j<colnum; ++j) { + for (i=0; i<rownum; ++i) { + y[j] += A[m*col_ids[j]+row_ids[i]]*x[i]; + } + } + + mxFree(col_ids); + } + + else { + + /* Do the actual computation */ + for (j=0; j<n; ++j) { + for (i=0; i<rownum; ++i) { + y[j] += A[m*j+row_ids[i]]*x[i]; + } + } + } + + + mxFree(row_ids); + + return; +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/rowlincomb.m Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,26 @@ +%ROWLINCOMB Linear combination of matrix rows. +% Y = ROWLINCOMB(X,A,ROWS) computes a linear combination of the rows of +% the matrix A. The row indices are specified in the vector ROWS, and the +% correspoinding coefficients are specified in the vector X. The vectors +% ROWS and X must be of the same length. The call Y = ROWLINCOMB(X,A,ROWS) +% is essentially equivalent to the command +% +% Y = X'*A(ROWS,:) . +% +% However, it is implemented much more efficiently. +% +% Y = ROWLINCOMB(X,A,ROWS,COLS) only works on the columns of A specified +% in COLS, returning a vector of length equal to COLS. This call is +% essentially equivalent to the command +% +% Y = X'*A(ROWS,COLS) . +% +% See also COLLINCOMB. + + +% Ron Rubinstein +% Computer Science Department +% Technion, Haifa 32000 Israel +% ronrubin@cs +% +% April 2009
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/sampgrid.m Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,72 @@ +function y = sampgrid(x,blocksize,varargin) +%SAMPGRID Sample a multi-dimensional matrix on a regular grid. +% Y = SAMPGRID(X,BLOCKSIZE,I1,I2,...,Ip) extracts block samples of size +% BLOCKSIZE from the p-dimensional matrix X, arranging the samples as the +% column vectors of the matrix Y. The locations of the (1,1,..,1)-th +% elements of each block are given in the index vectors I1,I2,..Ip. The +% total number of samples taken is length(I1)xlength(I2)x...xlength(Ip). +% BLOCKSIZE should either be a p-element vector of the form [N1,N2,...Np], +% or a scalar N which is shorthand for the square block size [N N ... N]. +% +% Example: Sample a set of blocks uniformly from a 2D image. +% +% n = 512; blocknum = 20000; blocksize = [8 8]; +% im = rand(n,n); +% [i1,i2] = reggrid(size(im)-blocksize+1, blocknum); +% blocks = sampgrid(im, blocksize, i1, i2); +% +% See also REGGRID. + +% Ron Rubinstein +% Computer Science Department +% Technion, Haifa 32000 Israel +% ronrubin@cs +% +% November 2007 + + +p = ndims(x); +if (p==2 && any(size(x)==1) && length(blocksize)==1) + p = 1; +end + +if (numel(blocksize)==1) + blocksize = ones(1,p)*blocksize; +end + +n = zeros(1,p); +for i = 1:p + n(i) = length(varargin{i}); +end + +nsamps = prod(n); + +% create y of the same class as x +y = zeros(prod(blocksize),nsamps,class(x)); + +% ids() contains the index of the current block in I1..Ip +ids = ones(p,1); + +% block_ids contains the indices of the current block in X +block_ids = cell(p,1); +for j = 1:p + block_ids{j} = varargin{j}(1) : varargin{j}(1)+blocksize(j)-1; +end + +for k = 1:nsamps + block = x(block_ids{:}); + y(:,k) = block(:); + + % increment ids() and block_ids{} + if (k<nsamps) + j = 1; + while (ids(j) == n(j)) + ids(j) = 1; + block_ids{j} = varargin{j}(1) : varargin{j}(1)+blocksize(j)-1; + j = j+1; + end + ids(j) = ids(j)+1; + block_ids{j} = varargin{j}(ids(j)) : varargin{j}(ids(j))+blocksize(j)-1; + end +end +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/scalarToRGB.m Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,18 @@ +function s = scalarToRGB(x,colors) +% input values are assumed to lie between 0 and 1 + +% Copyright 2008, Ewout van den Berg and Michael P. Friedlander +% http://www.cs.ubc.ca/labs/scl/sparco +% $Id: scalarToRGB.m 1040 2008-06-26 20:29:02Z ewout78 $ + +l = size(colors,1); +m = size(x,1); +n = size(x,2); +s = zeros(m,n,3); + +for i=1:m + for j=1:n + idx = max(1,min(l,1+floor((l-1) * x(i,j)))); + s(i,j,:) = colors(idx,:); + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/secs2hms.m Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,28 @@ +function [h,m,s] = secs2hms(t) +%SECS2HMS Convert seconds to hours, minutes and seconds. +% [H,M,S] = SECS2HMS(T) converts the specified number of seconds T to +% hours, minutes and seconds. H and M are whole numbers, and S is real. +% +% Example: Estimate the remaining time of a loop +% +% n = 10; tic; +% for i = 1:n +% pause(1); +% [h,m,s] = secs2hms( (n-i)*toc/i ); +% printf('estimated remaining time: %02d:%02d:%05.2f',h,m,s); +% end + + +% Ron Rubinstein +% Computer Science Department +% Technion, Haifa 32000 Israel +% ronrubin@cs +% +% April 2008 + + +s = t; +h = fix(s/3600); +s = rem(s,3600); +m = fix(s/60); +s = rem(s,60);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/spdiag.m Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,38 @@ +function Y = spdiag(V,K) +%SPDIAG Sparse diagonal matrices. +% SPDIAG(V,K) when V is a vector with N components is a sparse square +% matrix of order N+ABS(K) with the elements of V on the K-th diagonal. +% K = 0 is the main diagonal, K > 0 is above the main diagonal and K < 0 +% is below the main diagonal. +% +% SPDIAG(V) is the same as SPDIAG(V,0) and puts V on the main diagonal. +% +% See also DIAG, SPDIAGS. + + +% Ron Rubinstein +% Computer Science Department +% Technion, Haifa 32000 Israel +% ronrubin@cs +% +% June 2008 + + +if (nargin<2) + K = 0; +end + +n = length(V) + abs(K); + +if (K>0) + i = 1:length(V); + j = K+1:n; +elseif (K<0) + i = -K+1:n; + j = 1:length(V); +else + i = 1:n; + j = 1:n; +end + +Y = sparse(i,j,V(:),n,n); \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/sprow.c Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,94 @@ +/************************************************************************** + * + * File name: sprow.c + * + * Ron Rubinstein + * Computer Science Department + * Technion, Haifa 32000 Israel + * ronrubin@cs + * + * Last Updated: 24.8.2009 + * + *************************************************************************/ + + +#include "mex.h" +#include "mexutils.h" + + +/* Input Arguments */ + +#define A_IN prhs[0] +#define J_IN prhs[1] + + +/* Output Arguments */ + +#define X_OUT plhs[0] +#define ID_OUT plhs[1] + + +void mexFunction(int nlhs, mxArray *plhs[], + int nrhs, const mxArray*prhs[]) + +{ + double *pr, *x, *id, rowid; + mwIndex *ir, *jc; + mwSize m, n; + mwIndex i, j, k, l, rowlen; + + if (nrhs != 2) { + mexErrMsgTxt("GETSPROW requires two input arguments."); + } else if (nlhs > 2) { + mexErrMsgTxt("Too many output arguments."); + } + + checkmatrix(A_IN, "GETSPROW", "A"); + checksparse(A_IN, "GETSPROW", "A"); + checkscalar(J_IN, "GETSPROW", "J"); + + m = mxGetM(A_IN); + n = mxGetN(A_IN); + + rowid = mxGetScalar(J_IN); + if (rowid < 0) { + mexErrMsgTxt("Invalid row index."); + } + j = (mwIndex)(rowid + 1e-2); + if (j<1 || j>m) { + mexErrMsgTxt("Row index is out of range."); + } + j--; + + pr = mxGetPr(A_IN); + ir = mxGetIr(A_IN); + jc = mxGetJc(A_IN); + + /* Determine length of row */ + rowlen = 0; + for (i=0; i<jc[n]; ++i) { + rowlen += (ir[i]==j); + } + + /* Allocate output parameters */ + X_OUT = mxCreateDoubleMatrix(1, rowlen, mxREAL); + ID_OUT = mxCreateDoubleMatrix(1, rowlen, mxREAL); + + x = mxGetPr(X_OUT); + id = mxGetPr(ID_OUT); + + /* Compute j-th row */ + k=0; + for (l=1; l<=n; ++l) { + i = jc[l-1]; + while (i<jc[l] && ir[i]<j) { + i++; + } + if (i<jc[l] && ir[i]==j) { + x[k] = pr[i]; + id[k] = l; + k++; + } + } + +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/sprow.m Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,19 @@ +%SPROW Extract row of sparse matrix. +% X = SPROW(A,J) where A is a sparse matrix, returns the nonzero values in +% the row A(J,:). +% +% [X,ID] = SPROW(A,J) also returns the column indices of the nonzeros. +% +% Note that the call [X,ID] = SPROW(A,J) is equivalent (but more efficient +% than) the Matlab code +% +% IDS = find(A(J,:)); +% X = A(J,IDS); + + +% Ron Rubinstein +% Computer Science Department +% Technion, Haifa 32000 Israel +% ronrubin@cs +% +% August 2009
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/thumbFromOp.m Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,42 @@ +function P = thumbFromOp(op,m,n,sm,sn,grayscale) +% Output matrix P, of size m x n, sampled at +% sm x sn upper top + +% Copyright 2008, Ewout van den Berg and Michael P. Friedlander +% http://www.cs.ubc.ca/labs/scl/sparco +% $Id: thumbFromOp.m 1040 2008-06-26 20:29:02Z ewout78 $ + +if nargin < 6, grayscale = 0; end + +info = op([],0); + +sm = min(info{1},sm); +sn = min(info{2},sn); +M = zeros(sm,sn); + +for i=1:sn + v = zeros(info{2},1); v(i) = 1; + w = real(op(v,1)); + + M(:,i) = w(1:sm); +end + +mn = min(min(M)); +mx = max(max(M)); +M = (M - mn) / (mx-mn); + +idxm = floor(linspace(1,sm+1,m+1)); idxm = idxm(1:end-1); +idxn = floor(linspace(1,sn+1,n+1)); idxn = idxn(1:end-1); + +if grayscale + P = 1-M(idxm,idxn); +else + clrmap = hsv; + M = 1 + round(M * (length(clrmap)-1)); + P = zeros(m,n,3); + for j1=1:m + for j2=1:n + P(j1,j2,:) = clrmap(M(idxm(j1),idxn(j2)),:); + end + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/thumbPlot.m Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,38 @@ +function P = thumbPlot(P,x,y,color) + +% Copyright 2008, Ewout van den Berg and Michael P. Friedlander +% http://www.cs.ubc.ca/labs/scl/sparco +% $Id: thumbPlot.m 1040 2008-06-26 20:29:02Z ewout78 $ + +m = size(P,1); +n = size(P,2); +if (size(P,3) == 0) & (length(color) == 3) + % Convert to gray-scale + color = 0.30*color(1) + 0.59*color(2) + 0.11*color(3); +end + +mnx = min(x); % Minimum x +mxx = max(x); % Maximum x +mny = min(y); % Minimum y +mxy = max(y); % Maximum y +dy = (mxy - mny) * 0.1; % Offset on vertical axis +sx = (mxx - mnx) * 1.0; % Scale of horizontal axis +sy = (mxy - mny) * 1.2; % Scale of vertical axis + +if (sx < 1e-6), sx = 1; end +if (sy < 1e-6), sy = 1; end + +for i=1:length(x)-1 + x0 = floor(1 + (n-1) * (x(i ) - mnx) / sx); + x1 = floor(1 + (n-1) * (x(i+1) - mnx) / sx); + y0 = floor( (n-1) * (y(i ) - mny + dy) / sy); + y1 = floor( (n-1) * (y(i+1) - mny + dy) / sy); + + samples = 1+2*max(abs(x1-x0)+1,abs(y1-y0)+1); + c = linspace(0,1,samples); + idx = round((1-c)*x0 + c*x1); + idy = n - round((1-c)*y0 + c*y1); + for j=1:samples + P(idy(j),idx(j),:) = color; + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/thumbwrite.m Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,18 @@ +function thumbwrite(data,name,opts) + + +% Copyright 2008, Ewout van den Berg and Michael P. Friedlander +% http://www.cs.ubc.ca/labs/scl/sparco +% $Id: thumbwrite.m 1040 2008-06-26 20:29:02Z ewout78 $ + +% +% data Thumbnail data in range 0-1 +% name Name of file (no extension) +% opts +% .thumbtype Type of image (png, eps, ps, ...) +% .thumbdir Output directory +% + +[type,ext] = getFigureExt(opts.thumbtype); +data = round(data * 255) / 255; +imwrite(data,[opts.thumbpath,name,'.',ext],type);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/timerclear.m Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,37 @@ +function timerclear() +%TIMERCLEAR Clear all timers. +% TIMERCLEAR clears all currenly registered timers, invalidating all +% timer ids. +% +% Note: since registered timers do not consume CPU power except for when +% the TIMER<*> functions are called, this function is only useful in +% situations where a large number of timers have been initialized, and +% there is a need to reclaim memory. +% +% See also TIMERINIT, TIMERETA. + + +% Ron Rubinstein +% Computer Science Department +% Technion, Haifa 32000 Israel +% ronrubin@cs +% +% June 2008 + + +global utiltbx_timer_start_times % start times +global utiltbx_time_lastdisp % last display times +global utiltbx_timer_iternums % iteration numbers +global utiltbx_timer_lastiter % last queried iteration numbers +global utiltbx_timer_name % timer names +global utiltbx_timer_callfun % timer calling functions + + +% clear all timers % + +utiltbx_timer_start_times = []; +utiltbx_time_lastdisp = []; +utiltbx_timer_iternums = []; +utiltbx_timer_lastiter = []; +utiltbx_timer_name = []; +utiltbx_timer_callfun = [];
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/timereta.m Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,98 @@ +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 +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/timerinit.m Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,110 @@ +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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DL/RLS-DLA/private/updateFigure.m Tue Mar 15 12:20:59 2011 +0000 @@ -0,0 +1,93 @@ +function updateFigure(opts, figTitle, figFilename) + +% Copyright 2008, Ewout van den Berg and Michael P. Friedlander +% http://www.cs.ubc.ca/labs/scl/sparco +% $Id: updateFigure.m 1040 2008-06-26 20:29:02Z ewout78 $ + +% Ensure default values are available +opts.linewidth = getOption(opts,'linewidth', []); +opts.fontsize = getOption(opts,'fontsize', []); +opts.markersize = getOption(opts,'markersize',[]); + +% Output the plots +if opts.update + % Set the line width, font size and marker size + chld = [gca; get(gca,'Children')]; + lnwd = ones(length(chld),1) * NaN; + fnts = ones(length(chld),1) * NaN; + mrks = ones(length(chld),1) * NaN; + for i=1:length(chld) + conf = get(chld(i)); + if ~isempty(opts.linewidth) && isfield(conf,'LineWidth') + lnwd(i) = get(chld(i),'LineWidth'); + if (lnwd(i) == 0.5) % Default + set(chld(i),'Linewidth',opts.linewidth); + end + end + if ~isempty(opts.fontsize) && isfield(conf,'FontSize') + fnts(i) = get(chld(i),'FontSize'); + if (fnts(i) == 10) % Default + set(chld(i),'FontSize',opts.fontsize); + end + end + if ~isempty(opts.markersize) && isfield(conf,'MarkerSize') + mrks(i) = get(chld(i),'MarkerSize'); + if (mrks(i) == 6) % Default + set(chld(i),'MarkerSize',opts.markersize); + end + end + end + + for i=1:length(opts.figtype) + updateFigureType(opts.update, 0, opts.figtype{i}, ... + opts.figpath, figTitle, figFilename); + end + + % Restore the line-widths, font size + for i=1:length(chld) + if ~isnan(lnwd(i)) + set(chld(i),'LineWidth',lnwd(i)); + end + if ~isnan(fnts(i)) + set(chld(i),'FontSize',fnts(i)); + end + if ~isnan(mrks(i)) + set(chld(i),'MarkerSize',mrks(i)); + end + end + +end + +% Show the plot +if opts.show + updateFigureType(0,opts.show,'','',figTitle,''); +end + + + +function updateFigureType(update,show,figtype,figpath,figTitle,figFilename) +filename = [figpath,figFilename]; + +switch lower(figtype) + case {'pdf'} + cmdPostprocess = sprintf('!pdfcrop %s.pdf %s.pdf >& /dev/null', ... + filename, filename); + otherwise + cmdPostprocess = []; +end + +[figtype,figext] = getFigureExt(figtype); + +% Print the figure for output (without title) +if update + evalc(sprintf('print -d%s %s.%s;', figtype, filename, figext)); + + if ~isempty(cmdPostprocess) + eval(cmdPostprocess); + end +end + +% Add title if needed +if show + title(figTitle); +end