Mercurial > hg > smallbox
changeset 1:7750624e0c73 version0.5
(none)
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/README.txt Thu Nov 05 16:36:01 2009 +0000 @@ -0,0 +1,99 @@ + +SMALLbox Version 0.5 + +31st July, 2009 + + +--------------------------------------------------------------------------- + +Copyright (2009): Ivan Damnjanovic, Matthew Davies + +SMALLbox is distributed under the terms of the GNU General Public License 3 + +This program is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with this program. If not, see <http://www.gnu.org/licenses/>. + +--------------------------------------------------------------------------- + +SMALLbox is a MATLAB Evaluation Framework for the EU FET-OPEN Project, no: 225913: +Sparse Models, Algorithms, and Learning for Large-scale data (SMALL). + +The project team includes researchers from the following institutions: + +INRIA: http://www.irisa.fr/metiss/gribonval/ +University of Edinburgh: http://www.see.ed.ac.uk/~mdavies4/ +Queen Mary, University of London: http://www.elec.qmul.ac.uk/people/markp/ +EPFL: http://people.epfl.ch/pierre.vandergheynst +Technion: http://www.cs.technion.ac.il/~elad/ + +--------------------------------------------------------------------------- + +Version 0.5 of SMALLbox will download and install the following existing toolboxes +related to Sparse Representations, Compressed Sensing and Dictionary Learning: + + +Sparco: http://www.cs.ubc.ca/labs/scl/sparco/ + +SPGL1: http://www.cs.ubc.ca/labs/scl/spgl1/ + +SparseLab: http://sparselab.stanford.edu/ + +Sparsify: http://www.see.ed.ac.uk/~tblumens/sparsify/sparsify.html + +GPSR: http://www.lx.it.pt/~mtf/GPSR/ + +OMPbox and KSVDBox: http://www.cs.technion.ac.il/~ronrubin/ + +IMPORTANT: +In order to use SparseLab please register at +http://cgi.stanford.edu/group/sparselab/cgi-bin/register.pl + +IMPORTANT: +To successfully install all toolboxes you will need to have MEX setup to compile C files. +If this is not already setup, run "mex -setup" or type "help mex" in the MATLAB command prompt. + +IMPORTANT: +Because the toolboxes are downloaded automatically, you must have an internet connection +to successfully install SMALLbox. + +IMPORTANT: +If you are running Matlab on MAC OSX or Linux, you must start Matlab with the jvm enabled. +Not doing so, will prevent you being able to unzip the downloaded toolboxes. + +To install the toolbox run the command "SMALLboxsetup" from the MATLAB command prompt. + +Once installed, there are two optional demo functions that can be run, + +SMALL_solver_test.m : + +Example test of solvers from different toolboxes on Sparco compressed sensing problems + +and + +SMALL_solver_test_Audio.m : + +Example test of solvers on Sparco audio source separation problems + +--------------------------------------------------------------------------- + +For more information on the SMALL Project, please visit the following website: + +http://small.inria.fr + + +Contacts: ivan.damnjanovic@elec.qmul.ac.uk + matthew.davies@elec.qmul.ac.uk + + +This code is in experimental stage; any comments or bug reports are +very welcome.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/SMALLboxSetup.m Thu Nov 05 16:36:01 2009 +0000 @@ -0,0 +1,295 @@ +function SMALLboxSetup(varargin) +% +% SMALLboxSetup +% +% Will automatically download and install existing toolboxes +% on sparse representations and dictionary learning +% +% For this function an internet connection is required. +% +% SMALLbox initialisation +% Ivan Damnjanovic, Matthew Davies 2009 + + +clc; + +FS=filesep; + +fprintf('\n ********************************************************************'); +fprintf('\n\n This script will install the SMALLbox Evaluation Framework v.0.4'); +fprintf('\n\n It contains the following toolboxes:'); +fprintf('\n Sparco version 1.2, incorporating Rice Wavelet Toolbox version 2.4'); +fprintf('\n SPGL1 Toolbox version 1.7 '); +fprintf('\n SparseLab Toolbox version 2.1'); +fprintf('\n Sparsify Toolbox version 0.4'); +fprintf('\n GPSR Toolbox version 5.0'); +fprintf('\n OMPbox version 9'); +fprintf('\n KSVDbox version 10'); +fprintf('\n\n ********************************************************************'); + +fprintf('\n\n The toolbox will be installed in: '); +fprintf('\n %s%s\n',pwd,FS); +fprintf('\n ********************************************************************'); +fprintf('\n\n IMPORTANT: To successfully install all toolboxes'); +fprintf('\n you will need to have MEX setup to compile C files.'); +fprintf('\n\n If this is not already setup, please type "n" to exit and then '); +fprintf('\n run "mex -setup" or type "help mex" in the MATLAB command prompt.'); +fprintf('\n\n ********************************************************************'); + + +fprintf('\n ********************************************************************'); +fprintf('\n\n IMPORTANT: YOU MUST HAVE AN INTERNET CONNECTION'); +fprintf('\n YOU CANNOT INSTALL SMALLBOX WITHOUT ONE!'); +fprintf('\n\n ********************************************************************'); +install_ack = input('\n\n Do you wish to continue: ([y]/n)? ','s'); + +if strcmp(install_ack,'"n"'), + install_ack = 'n'; +end + +if install_ack == 'n', + return; +else + fprintf('\n\n Installation now beginning...'); +% fprintf('\n\nFor further information on the toolboxes see SMALLbox documentation\n\n'); +end + + + + +SMALL_path=pwd; +SMALL_p=genpath(SMALL_path); +addpath(SMALL_p); + + +fprintf('\n ******************************************************************'); +fprintf('\n\n Initialising SPARCO and Rice Wavelet Toolbox Setup'); + + + Sparco_path = [SMALL_path,FS,'toolboxes',FS,'SPARCO']; + if exist([Sparco_path, FS, 'sparco-1.2.zip'],'file'), + Sparco_zip=[Sparco_path, FS, 'sparco-1.2.zip']; + else + Sparco_zip='http://www.cs.ubc.ca/labs/scl/sparco/downloads.php?filename=sparco-1.2.zip'; + fprintf('\n\n Downloading toolbox, please be patient\n\n'); + end + unzip(Sparco_zip,Sparco_path); + Sparco_p=genpath(Sparco_path); + addpath(Sparco_p); + cd(SMALL_path) + + + + if exist('curvelab.pdf','file') + crvroot = fileparts(which('curvelab.pdf')); + addtopath(crvroot,'fdct_usfft_matlab'); + addtopath(crvroot,'fdct_wrapping_matlab'); + addtopath(crvroot,'fdct_wrapping_cpp/mex'); + addtopath(crvroot,'fdct3d/mex'); + else + fprintf(['\nWarning: CurveLab is not in the path. Sparco Problems 50-51 ' ... + 'will not work.\n\n']); + end + + cd([Sparco_path, FS, 'sparco-1.2', FS, 'tools' ,FS, 'rwt']) + fprintf('Compiling the Rice Wavelet Toolbox MEX interfaces...'); + try + if exist('mdwt' ,'file')~=3, mex mdwt.c mdwt_r.c; end + if exist('midwt' ,'file')~=3, mex midwt.c midwt_r.c; end + if exist('mrdwt' ,'file')~=3, mex mrdwt.c mrdwt_r.c; end + if exist('mirdwt','file')~=3, mex mirdwt.c mirdwt_r.c; end + fprintf('SPARCO Installation Successful!\n'); + catch + warning('Could not compile Rice Wavelet Toolbox MEX interfaces.'); + end + cd(SMALL_path) + + + +fprintf('\n ******************************************************************'); +fprintf('\n\n Initialising SPGL1 Setup'); + + try + SPGL1_path = [SMALL_path,FS,'toolboxes',FS,'SPGL1']; + if exist([SPGL1_path, FS, 'spgl1-1.7.zip'],'file'), + SPGL1_zip=[SPGL1_path, FS, 'spgl1-1.7.zip']; + else + SPGL1_zip='http://www.cs.ubc.ca/labs/scl/spgl1/downloads.php?filename=spgl1-1.7.zip'; + fprintf('\n\n Downloading toolbox, please be patient\n\n'); + end + unzip(SPGL1_zip,SPGL1_path); + SPGL1_p=genpath(SPGL1_path); + addpath(SPGL1_p); + + + cd([SPGL1_path,FS,'spgl1-1.7']); + fprintf('Compiling SPGL1 MEX interfaces ...'); + try + spgsetup; + fprintf('\n SPGL1 Installation Successful!\n'); + catch + warning('Could not compile SPGL1 MEX interfaces.'); + end + catch + fprintf('\n SPGL1 Installation Failed\n'); + end + cd(SMALL_path); + +fprintf('\n ******************************************************************'); +fprintf('\n\n Initialising SparseLab Setup'); + + try + SL_path = [pwd,FS,'toolboxes',FS,'SparseLab']; + if exist([SL_path, FS, 'SparseLab21-Core.zip'],'file'), + SL_zip=[SL_path, FS, 'SparseLab21-Core.zip']; + else + SL_zip='http://sparselab.stanford.edu/SparseLab_files/Download_files/SparseLab21-Core.zip'; + fprintf('\n\n Downloading toolbox, please be patient\n\n'); + end + unzip(SL_zip,SL_path); + SL_p=genpath(SL_path); + addpath(SL_p); + fprintf('\n SparseLab Installation Successful!\n'); + catch + fprintf('\n SparseLab Installation Failed\n'); + cd(SMALL_path); + end + + +fprintf('\n ******************************************************************'); +fprintf('\n\n Initialising Sparsify Setup'); + + try + Sparsify_path = [pwd,FS,'toolboxes',FS,'Sparsify']; + if exist([Sparsify_path, FS, 'sparsify_0_4.zip'],'file'), + Sparsify_zip=[Sparsify_path, FS, 'sparsify_0_4.zip']; + else + Sparsify_zip='http://www.see.ed.ac.uk/~tblumens/sparsify/sparsify_0_4.zip'; + fprintf('\n\n Downloading toolbox, please be patient\n\n'); + end + unzip(Sparsify_zip,Sparsify_path); + Sparsify_p=genpath(Sparsify_path); + addpath(Sparsify_p); + fprintf('\n Sparsify Installation Successful\n'); + catch + fprintf('\n Sparsify Installation Failed\n'); + cd(SMALL_path); + end + +fprintf('\n ******************************************************************'); +fprintf('\n\n Initialising GPSR Setup'); + + try + GPSR_path = [pwd,FS,'toolboxes',FS,'GPSR']; + if exist([GPSR_path, FS, 'GPSR_6.0.zip'],'file'), + GPSR_zip=[GPSR_path, FS,'GPSR_6.0.zip']; + else + GPSR_zip='http://www.lx.it.pt/~mtf/GPSR/GPSR_6.0.zip'; + fprintf('\n\n Downloading toolbox, please be patient\n\n'); + end + unzip(GPSR_zip,GPSR_path); + GPSR_p=genpath(GPSR_path); + addpath(GPSR_p); + fprintf('\n GPSR Installation Successful\n'); + catch + fprintf('\n GPSR Installation Failed'); + cd(SMALL_path); + end + + +fprintf('\n ******************************************************************'); +fprintf('\n\n Initialising OMPbox and KSVDBox Setup'); + + try + KSVD_path = [pwd,FS,'toolboxes',FS,'KSVD']; + if exist([KSVD_path, FS, 'ompbox9.zip'],'file'), + omp_zip=[KSVD_path, FS, 'ompbox9.zip']; + else + omp_zip='http://www.cs.technion.ac.il/%7Eronrubin/Software/ompbox9.zip'; + fprintf('\n\n Downloading toolbox, please be patient\n\n'); + end + unzip(omp_zip,[KSVD_path, FS, 'ompbox']); + + cd([KSVD_path, FS, 'ompbox', FS, 'private']); + make; + cd(SMALL_path); + + if exist([KSVD_path, FS, 'ksvdbox10.zip'],'file'), + KSVD_zip=[KSVD_path, FS, 'ksvdbox10.zip']; + else + KSVD_zip='http://www.cs.technion.ac.il/%7Eronrubin/Software/ksvdbox10.zip'; + fprintf('\n\n Downloading toolbox, please be patient\n\n'); + end + unzip(KSVD_zip,[KSVD_path, FS, 'ksvdbox']); + cd([KSVD_path, FS, 'ksvdbox', FS, 'private']); + make; + cd(SMALL_path); + KSVD_p=genpath(KSVD_path); + addpath(KSVD_p); + fprintf('\n KSVDBox and OMPBox Installation Successful\n'); + catch + fprintf('\n KSVDBox and OMPBox Installation Failed'); + cd(SMALL_path); + end + + + + +fprintf('\n ******************************************************************'); +fprintf('\n\n Initialising SMALLbox Examples Setup'); + + % Need to do a bit of temporary housekeeping first. + cd(SMALL_path); + try + cd(['examples',FS,'private']); + if exist('addtocols' ,'file')~=3, + fprintf('\n Compiling MEX interfaces for SMALL examples \n'); + make; + end + fprintf('\n SMALLbox Examples Installation Successful! \n'); + catch + fprintf('\n SMALLbox Examples Installation Failed \n'); + end + cd(SMALL_path); + + + +fprintf('\n ******************************************************************'); +fprintf('\n\n SMALLbox Installation Complete!'); + + +fprintf('\n\n For more information on the installed toolboxes see'); +fprintf('\n\n Sparco: http://www.cs.ubc.ca/labs/scl/sparco/'); +fprintf('\n\n SPGL1: http://www.cs.ubc.ca/labs/spgl1/?n=HomePage'); +fprintf('\n\n SparseLab: http://sparselab.stanford.edu/ (PLEASE REGISTER SPARSELAB!)'); +fprintf('\n\n Sparsify: http://www.see.ed.ac.uk/~tblumens/sparsify/sparsify.html'); +fprintf('\n\n GPSR: http://www.lx.it.pt/~mtf/GPSR/'); +fprintf('\n\n OMPbox and KSVDBox: http://www.cs.technion.ac.il/~ronrubin/\n'); + + +% LIST DEMOS FROM EXAMPLE DIRECTORY... + +demo_ack = input('\n\n Would you like to run a demo: ([y]/n)? ','s'); + +if demo_ack == 'n', + fprintf('\n Thank you for installing SMALLbox.'); + fprintf('\n For information on the SMALLbox example scripts'); + fprintf('\n Please see the examples directory. \n'); + return; +else + + demo_choice = input('\n Enter 1 to run SMALL_solver_test, 2 to run SMALL_solver_test_Audio or q to quit: ','s'); + switch(demo_choice) + case{'1'} + fprintf('\n Running SMALL_solver_test problem'); + SMALL_solver_test; + case{'2'} + fprintf('\n Running SMALL_solver_test_Audio problem'); + SMALL_solver_test_Audio; + otherwise + return; + end + +end + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/examples/SMALL_solver_test.m Thu Nov 05 16:36:01 2009 +0000 @@ -0,0 +1,129 @@ +function SMALL_solver_test +% Example test of solvers from different toolboxes on Sparco compressed +% sensing problems +% +% The main purpose of this example is to show how to use SMALL structure +% to solve SPARCO compressed sensing problems (1-11) and compare results +% from different solvers. +% To generate SMALL.Problem part of structure you can use generateProblem +% function from Sparco toolbox giving the problem number and any +% additional parameters you might want to change. Alternatively, you can +% might want to consult sparco documentation to write a problem by +% yourself. There are four fields the must be specified in SMALL.Problem +% - A, b, sizeA and reconstruct. +% +% To generate SMALL.solver part of the structure you must specify three +% fields: +% +% SMALL.solver.toolbox - string with toolbox name is needed because +% different toolboxes are calling solver +% functions in different ways. +% SMALL.solver.name - its string representing solver name (e.g. +% SolveOMP) +% SMALL.solver.param - string that contains optional parameters for +% particular solver (all parameters you want to +% specify except A, b and size of solution) +% +% Every call to SMALL_solve function will generate following output: +% +% SMALL.solver.solution - contains solution vector x +% SMALL.solver.reconstructed - vector containing signal reconstructed +% from the solution +% SMALL.solver.time - time that solver spent to find the solution +% +% SMALL_plot function plots the SMALL.solver.solution and reconstructed +% against original signal. +% +% In this particular example we are testing SMALL_cgp, SMALL_chol, +% SolveOMP form SparseLab and greed_pcgp form Sparsify against "PROB006 +% Daubechies basis, Gaussian ensemble measurement basis, piecewise cubic +% polynomial signal" from Sparco. +% +% +% Ivan Damnjanovic 2009% +% +% SPARCO Copyright 2008, Ewout van den Berg and Michael P. Friedlander +% http://www.cs.ubc.ca/labs/scl/sparco +% $Id: exGPSR.m 1040 2008-06-26 20:29:02Z ewout78 $ + +fprintf('\n\nExample test of SMALL solver against their counterparts on Sparco problems.\n\n'); + +%% +% Generate SPARCO problem + +global SMALL +SMALL.Problem = generateProblem(6, 'P', 6, 'm', 270,'n',1024, 'show'); +%% + +%% +% SMALL Conjugate Gradient test + +SMALL.solver.toolbox='SMALL'; +SMALL.solver.name='SMALL_cgp'; + +% In the following string all parameters except matrix, measurement vector +% and size of solution need to be specified. If you are not sure which +% parameters are needed for particular solver type "help <Solver name>" in +% MATLAB command line + +SMALL.solver.param='200, 1e-14'; + +SMALL=SMALL_solve(SMALL); + + +SMALL_plot(SMALL); + +%% +% SMALL OMP with Cholesky update test + +SMALL.solver.toolbox='SMALL'; +SMALL.solver.name='SMALL_chol'; + +% In the following string all parameters except matrix, measurement vector +% and size of solution need to be specified. If you are not sure which +% parameters are needed for particular solver type "help <Solver name>" in +% MATLAB command line + +SMALL.solver.param='200, 1e-14'; + +SMALL=SMALL_solve(SMALL); + +SMALL_plot(SMALL); +%% +% SolveOMP from SparseLab test + +SMALL.solver.toolbox='SparseLab'; +SMALL.solver.name='SolveOMP'; + +% In the following string all parameters except matrix, measurement vector +% and size of solution need to be specified. If you are not sure which +% parameters are needed for particular solver type "help <Solver name>" in +% MATLAB command line + +SMALL.solver.param='200, 0, 0, 0, 1e-14'; + +SMALL=SMALL_solve(SMALL); + +SMALL_plot(SMALL); + +%% +% greed_pcgp from Sparsify test + +SMALL.solver.toolbox='Sparsify'; +SMALL.solver.name='greed_pcgp'; + +% In the following string all parameters except matrix, measurement vector +% and size of solution need to be specified. If you are not sure which +% parameters are needed for particular solver type "help <Solver name>" in +% MATLAB command line + +SMALL.solver.param='''stopCrit'', ''M'', ''stopTol'', 200'; + +SMALL=SMALL_solve(SMALL); + +SMALL_plot(SMALL); + +%% + + +end % function SMALL_solver_test
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/examples/SMALL_solver_test_Audio.m Thu Nov 05 16:36:01 2009 +0000 @@ -0,0 +1,93 @@ +function SMALL_solver_test_Audio +% Example test of solvers on Sparco audio source separation problems +% +% The main purpose of this example is to show how to use SMALL structure +% to solve SPARCO audio source3 separation problems (401-402) and to +% compare results from different solvers. +% To generate SMALL.Problem part of structure you can use generateProblem +% function from Sparco toolbox giving the problem number and any +% additional parameters you might want to change. Alternatively, you can +% might want to consult sparco documentation to write a problem by +% yourself. There are four fields the must be specified in SMALL.Problem +% - A, b, sizeA and reconstruct. +% +% To generate SMALL.solver part of the structure you must specify three +% fields: +% +% SMALL.solver.toolbox - string with toolbox name is needed because +% different toolboxes are calling solver +% functions in different ways. +% SMALL.solver.name - its string representing solver name (e.g. +% SolveBP) +% SMALL.solver.param - string that contains optional parameters for +% particular solver (all parameters you want to +% specify except A, b and size of solution) +% +% Every call to SMALL_solve function will generate following output: +% +% SMALL.solver.solution - contains solution vector x +% SMALL.solver.reconstructed - vector containing signal reconstructed +% from the solution +% SMALL.solver.time - time that solver spent to find the solution +% +% SMALL_plot function plots the SMALL.solver.solution and reconstructed +% sources against original audio sources. +% SMALL_playAudio function plays audio sources of original and +% reconstructed signal as well as mixed signal. +% +% +% Ivan Damnjanovic 2009% +% +% SPARCO Copyright 2008, Ewout van den Berg and Michael P. Friedlander +% http://www.cs.ubc.ca/labs/scl/sparco +% $Id: exGPSR.m 1040 2008-06-26 20:29:02Z ewout78 $ + +fprintf('\n\nExample test of solvers on Sparco Audio problems (401,402).\n\n'); + +%% +% Generate SPARCO problem + +global SMALL +SMALL.Problem = generateProblem(402,'show'); +%% + +%% +% SMALL Conjugate Gradient test + +SMALL.solver.toolbox='SMALL'; +SMALL.solver.name='SMALL_cgp'; + +% In the following string all parameters except matrix, measurement vector +% and size of solution need to be specified. If you are not sure which +% parameters are needed for particular solver type "help <Solver name>" in +% MATLAB command line + +SMALL.solver.param='1500, 1e-14'; + +SMALL=SMALL_solve(SMALL); + + +SMALL_plot(SMALL); +SMALL_playAudio(SMALL); + +%% +% SolveOMP from SparseLab test + +SMALL.solver.toolbox='SparseLab'; +SMALL.solver.name='SolveBP'; + +% In the following string all parameters except matrix, measurement vector +% and size of solution need to be specified. If you are not sure which +% parameters are needed for particular solver type "help <Solver name>" in +% MATLAB command line + +SMALL.solver.param='10'; + +SMALL=SMALL_solve(SMALL); + +SMALL_plot(SMALL); +SMALL_playAudio(SMALL); +%% + + +end % function SMALL_solver_test
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/examples/private/add_dc.m Thu Nov 05 16:36:01 2009 +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/examples/private/addtocols.c Thu Nov 05 16:36:01 2009 +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/examples/private/addtocols.m Thu Nov 05 16:36:01 2009 +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/examples/private/collincomb.c Thu Nov 05 16:36:01 2009 +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/examples/private/collincomb.m Thu Nov 05 16:36:01 2009 +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/examples/private/countcover.m Thu Nov 05 16:36:01 2009 +0000 @@ -0,0 +1,31 @@ +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. + + +% 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/examples/private/diag_ids.m Thu Nov 05 16:36:01 2009 +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/examples/private/dictdist.m Thu Nov 05 16:36:01 2009 +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/examples/private/imnormalize.m Thu Nov 05 16:36:01 2009 +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/examples/private/iswhole.m Thu Nov 05 16:36:01 2009 +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/examples/private/make.m Thu Nov 05 16:36:01 2009 +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 +% +% April 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'}; + +for i = 1:length(sourcefiles) + printf('Compiling %s...', sourcefiles{i}); + mex(sourcefiles{i},compile_params{:}); +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/examples/private/normcols.m Thu Nov 05 16:36:01 2009 +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/examples/private/printf.m Thu Nov 05 16:36:01 2009 +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/examples/private/reggrid.m Thu Nov 05 16:36:01 2009 +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/examples/private/remove_dc.m Thu Nov 05 16:36:01 2009 +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/examples/private/rowlincomb.c Thu Nov 05 16:36:01 2009 +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/examples/private/rowlincomb.m Thu Nov 05 16:36:01 2009 +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/examples/private/sampgrid.m Thu Nov 05 16:36:01 2009 +0000 @@ -0,0 +1,93 @@ +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 (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 + + +% +% p = ndims(x); +% +% 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)); +% +% id = cell(p,1); +% for k = 1:nsamps +% [id{:}] = ind2sub(n,k); +% for j = 1:p +% id{j} = varargin{j}(id{j}) : varargin{j}(id{j})+blocksize(j)-1; +% end +% block = x(id{:}); +% y(:,k) = block(:); +% end +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/examples/private/secs2hms.m Thu Nov 05 16:36:01 2009 +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/examples/private/spdiag.m Thu Nov 05 16:36:01 2009 +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/examples/private/timerclear.m Thu Nov 05 16:36:01 2009 +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/examples/private/timereta.m Thu Nov 05 16:36:01 2009 +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/examples/private/timerinit.m Thu Nov 05 16:36:01 2009 +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/solvers/SMALL_MP.m Thu Nov 05 16:36:01 2009 +0000 @@ -0,0 +1,85 @@ +function [A]=SMALL_MP(Dict,X, m, maxNumCoef, errorGoal, varargin) +%% +%============================================= +% Sparse coding of a group of signals based on a given +% dictionary and specified number of atoms to use. +% input arguments: Dict - the dictionary +% X - the signals to represent +% m - number of atoms in Dictionary +% errorGoal - the maximal allowed representation error for +% each signal. +% +% optional: if Dict is function handle then Transpose Dictionary +% handle needs to be specified. +% +% based on KSVD toolbox solver found on Miki Elad webpage (finding inverse +% with pinv() is changed with matching pursuit) +% Ivan Damnjanovic 2009 +%============================================= +%% +%% +% This Dictionary check is based on Thomas Blumensath work in sparsify 0_4 greedy solvers +explicitD=0; +if isa(Dict,'float') + D =@(z) Dict*z; + Dt =@(z) Dict'*z; + explicitD=1; +elseif isobject(Dict) + D =@(z) Dict*z; + Dt =@(z) Dict'*z; +elseif isa(Dict,'function_handle') + try + DictT=varargin{1}; + if isa(DictT,'function_handle'); + D=Dict; + Dt=DictT; + else + error('If Dictionary is a function handle,Transpose Dictionary also needs to be a function handle. '); + end + catch + error('If Dictionary is a function handle, Transpose Dictionary needs to be specified. Exiting.'); + end +else + error('Dictionary is of unsupported type. Use explicit matrix, function_handle or object. Exiting.'); +end +%% +[n,P]=size(X); + +E2 = errorGoal^2; + + +A = sparse(m,size(X,2)); +for k=1:1:P, + + x = X(:,k); + residual=x; + indx = []; + j=0; + + currResNorm = norm(residual); + errorGoal=errorGoal*currResNorm; + + a = zeros(m,1); + + + while currResNorm>errorGoal & j < maxNumCoef, + + j = j+1; + dir=Dt(residual); + + [tmp__, pos]=max(abs(dir)); + indx(j)=pos; + a(pos)=dir(pos); + e = zeros(m,1); + e(pos) = 1; + residual=residual-D(e)*a(pos); + + currResNorm = norm(residual); + end; + if (~isempty(indx)) + A(indx,k)=a(indx); + end +end; +return; + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/solvers/SMALL_cgp.m Thu Nov 05 16:36:01 2009 +0000 @@ -0,0 +1,103 @@ +function [A]=SMALL_cgp(Dict,X, m, maxNumCoef, errorGoal, varargin) +%============================================= +% Sparse coding of a group of signals based on a given +% dictionary and specified number of atoms to use. +% input arguments: Dict - the dictionary +% X - the signals to represent +% m - number of atoms in Dictionary +% errorGoal - the maximal allowed representation error for +% each signal. +% +% optional: if Dict is function handle then Transpose Dictionary +% handle needs to be specified. +% +% output arguments: A - sparse coefficient matrix. +% +% based on KSVD toolbox solver found on Miki Elad webpage (finding inverse +% with pinv() is changed with conjugate gradient method) +% Ivan Damnjanovic 2009 +%============================================= +%% +%% + +% This Dictionary check is based on Thomas Blumensath work in sparsify 0_4 greedy solvers +if isa(Dict,'float') + D =@(z) Dict*z; + Dt =@(z) Dict'*z; +elseif isobject(Dict) + D =@(z) Dict*z; + Dt =@(z) Dict'*z; +elseif isa(Dict,'function_handle') + try + DictT=varargin{1}; + if isa(DictT,'function_handle'); + D=Dict; + Dt=DictT; + else + error('If Dictionary is a function handle,Transpose Dictionary also needs to be a function handle. '); + end + catch + error('If Dictionary is a function handle, Transpose Dictionary needs to be specified. Exiting.'); + end +else + error('Dictionary is of unsupported type. Use explicit matrix, function_handle or object. Exiting.'); +end +%% + +[n,P]=size(X); + + + +A = sparse(m,size(X,2)); + +for k=1:1:P, + + x = X(:,k); + residual = x; + indx =[]; + j=0; + + + currResNorm = residual'*residual/n; + errorGoal=errorGoal*currResNorm; + a = zeros(m,1); + p = zeros(m,1); + + while j < maxNumCoef, + + j = j+1; + dir=Dt(residual); + + [tmp__, pos]=max(abs(dir)); + indx(j)=pos; + + p(indx)=dir(indx); + Dp=D(p); + pDDp=Dp'*Dp; + + if (j==1) + beta=0; + else + beta=-Dp'*residual/pDDp; + end + + alfa=residual'*Dp/pDDp; + a=a+alfa*p; + p(indx)=dir(indx)+beta*p(indx); + + residual=residual-alfa*Dp; + + currResNorm = residual'*residual/n; + if currResNorm<errorGoal + fprintf('\nFound exact representation! \n'); + break; + end + end; + if (~isempty(indx)) + A(indx,k)=a(indx); + end +end; +return; + + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/solvers/SMALL_chol.m Thu Nov 05 16:36:01 2009 +0000 @@ -0,0 +1,130 @@ +function [A]=SMALL_chol(Dict,X, m, maxNumCoef, errorGoal, varargin) +%% +%============================================= +% Sparse coding of a group of signals based on a given +% dictionary and specified number of atoms to use. +% input arguments: Dict - the dictionary +% X - the signals to represent +% m - number of atoms in Dictionary +% errorGoal - the maximal allowed representation error for +% each signal. +% +% optional: if Dict is function handle then Transpose Dictionary +% handle needs to be specified. +% +% output arguments: A - sparse coefficient matrix. +% +% based on KSVD toolbox solver found on Miki Elad webpage (finding inverse +% with pinv() is changed with OMP Cholesky update) +% Ivan Damnjanovic 2009 +%============================================= +%% +% This Dictionary check is based on Thomas Blumensath work in sparsify 0_4 greedy solvers +explicitD=0; +if isa(Dict,'float') + D =@(z) Dict*z; + Dt =@(z) Dict'*z; + explicitD=1; +elseif isobject(Dict) + D =@(z) Dict*z; + Dt =@(z) Dict'*z; +elseif isa(Dict,'function_handle') + try + DictT=varargin{1}; + if isa(DictT,'function_handle'); + D=Dict; + Dt=DictT; + else + error('If Dictionary is a function handle,Transpose Dictionary also needs to be a function handle. '); + end + catch + error('If Dictionary is a function handle, Transpose Dictionary needs to be specified. Exiting.'); + end +else + error('Dictionary is of unsupported type. Use explicit matrix, function_handle or object. Exiting.'); +end +%% +[n,P]=size(X); + + + +global opts opts_tr machPrec +opts.UT = true; +opts_tr.UT = true; opts_tr.TRANSA = true; +machPrec = 1e-5; + +A = sparse(m,size(X,2)); +for k=1:1:P, + + R_I = []; + x=X(:,k); + residual=x; + indx = []; + a = zeros(m,1); + currResNorm = norm(residual); + errorGoal=errorGoal*currResNorm; + j = 0; + while currResNorm>errorGoal & j < maxNumCoef, + j = j+1; + dir=Dt(residual); + + [tmp__, pos]=max(abs(dir)); + + [R_I, flag] = updateChol(R_I, n, m, D, explicitD, indx, pos, Dt); + + + indx(j)=pos; + dx=zeros(m,1); + + z = linsolve(R_I,dir(indx),opts_tr); + + dx(indx) = linsolve(R_I,z,opts); + a(indx) = a(indx) + dx(indx); + + residual=x-D(a); + currResNorm = norm(residual); + + + end; + if (~isempty(indx)) + A(indx,k)=a(indx); + end +end; +return; + + +function [R, flag] = updateChol(R, n, N, A, explicitA, activeSet, newIndex, varargin) + +% updateChol: Updates the Cholesky factor R of the matrix +% A(:,activeSet)'*A(:,activeSet) by adding A(:,newIndex) +% If the candidate column is in the span of the existing +% active set, R is not updated, and flag is set to 1. + +global opts_tr machPrec +flag = 0; + +if (explicitA) + newVec = A(:,newIndex); +else + At=varargin{1}; + e = zeros(N,1); + e(newIndex) = 1; + newVec = A(e);%feval(A,1,n,N,e,1:N,N); +end + +if isempty(activeSet), + R = sqrt(sum(newVec.^2)); +else + if (explicitA) + p = linsolve(R,A(:,activeSet)'*A(:,newIndex),opts_tr); + else + AnewVec = At(newVec);%feval(A,2,n,length(activeSet),newVec,activeSet,N); + p = linsolve(R,AnewVec(activeSet),opts_tr); + end + q = sum(newVec.^2) - sum(p.^2); + if (q <= machPrec) % Collinear vector + flag = 1; + else + R = [R p; zeros(1, size(R,2)) sqrt(q)]; + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/util/SL_A.m Thu Nov 05 16:36:01 2009 +0000 @@ -0,0 +1,20 @@ + function y = SL_A(mode, m, n, x, I, dim) + % Ivan Damnjanovic 2009 + % This is auxilary function to allow implicit matrices from SPARCO + % to be used with SparsLab solvers + + global SMALL + if (mode == 1) + + u = zeros(dim, 1); + u(I) = x; + y = SMALL.Problem.A(u,1); + + elseif (mode == 2) + + x2 = SMALL.Problem.A(x,2); + y = x2(I); + + end + + end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/util/SMALL_playAudio.m Thu Nov 05 16:36:01 2009 +0000 @@ -0,0 +1,40 @@ +function SMALL_playAudio(SMALL) + % Ivan Damnjanovic 2009 + % Function gets as input SMALL structure and play the original and + % the reconstructed signal + + +ch=''; +while 1 + request = input('\nWhat do you want to hear? \n 1. Original signal \n 2. Mixed \n 3. Reconstructed signal \n 4. Quit player\n','s'); + request = sscanf(request,'%d'); + switch request + case 1 + fprintf('Original signal has %d sources.\n Which one do you want to hear?', size(SMALL.Problem.signal,2)) + fprintf('Enter a number between 1 and %d', size(SMALL.Problem.signal,2)) + ch=input('','s'); + ch=sscanf(ch,'%d'); + if (ch>=1)&&(ch<=size(SMALL.Problem.signal,2)) + soundsc(SMALL.Problem.signal(:,ch),8000); + end + case 2 + fprintf('Number of channels in mixed signal is %d.\n Which one do you want to hear?', size(SMALL.Problem.b,2)) + fprintf('Enter a number between 1 and %d', size(SMALL.Problem.signal,2)) + ch=input('','s'); + ch=sscanf(ch,'%d'); + if (ch>=1)&&(ch<=size(SMALL.Problem.b,2)) + soundsc(SMALL.Problem.b(:,ch),8000); + end + case 3 + fprintf('Reconstructed signal has %d sources.\n Which one do you want to hear?', size(SMALL.solver.reconstructed,2)) + fprintf('Enter a number between 1 and %d', size(SMALL.Problem.signal,2)) + ch=input('','s'); + ch=sscanf(ch,'%d'); + if (ch>=1)&&(ch<=size(SMALL.solver.reconstructed,2)) + soundsc(SMALL.solver.reconstructed(:,ch),8000); + end + case 4 + return; + end +end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/util/SMALL_plot.m Thu Nov 05 16:36:01 2009 +0000 @@ -0,0 +1,23 @@ +function SMALL_plot(SMALL) + % Ivan Damnjanovic 2009 + % Function gets as input SMALL structure and plots the solution and + % reconstructed signal + + figure; + + plot(1:length(SMALL.solver.solution), SMALL.solver.solution, 'b'); + title(['Coefficients of the reconstructed signal using ', SMALL.solver.name, ' with parameters ', SMALL.solver.param]) + xlabel('Coefficient') + + % Use SMALL.solution to reconstruct the signal. + + % Plot reconstructed signal against original + + + for i=1:size(SMALL.solver.reconstructed,2) + figure; + plot(1:length(SMALL.solver.reconstructed(:,i)), SMALL.solver.reconstructed(:,i) ,'b.-', 1:length(SMALL.Problem.signal(:,i)), SMALL.Problem.signal(:,i),'r--'); + legend(SMALL.solver.name,'Original signal'); + title('Reconstructed and original signals'); + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/util/SMALL_solve.m Thu Nov 05 16:36:01 2009 +0000 @@ -0,0 +1,32 @@ +function SMALL = SMALL_solve(SMALL) + % Ivan Damnjanovic 2009 + % Function gets as input SMALL structure that contains SPARCO problem to + % be solved, name of the toolbox and solver, and parameters file for particular solver + % Outputs are solution, reconstructed signal and time spent + %% + + A = @(x) SMALL.Problem.A(x,1); % The operator + AT = @(y) SMALL.Problem.A(y,2); % and its transpose. + b = SMALL.Problem.b; % The right-hand-side vector. + m = SMALL.Problem.sizeA(1); % m is the no. of rows. + n = SMALL.Problem.sizeA(2); % n is the no. of columns. + + fprintf('\nStarting solver %s... \n', SMALL.solver.name); + start=cputime; + %% + if strcmpi(SMALL.solver.toolbox,'sparselab') + y = eval([SMALL.solver.name,'(''SL_A'', b, n,',SMALL.solver.param,');']); + elseif strcmpi(SMALL.solver.toolbox,'sparsify') + y = eval([SMALL.solver.name,'(b, A, n, ''P_trans'', AT,',SMALL.solver.param,');']); + elseif (strcmpi(SMALL.solver.toolbox,'spgl1')||strcmpi(SMALL.solver.toolbox,'gpsr')) + y = eval([SMALL.solver.name,'(b, A,',SMALL.solver.param,');']); + else + y = eval([SMALL.solver.name,'(A, b, n,',SMALL.solver.param,',AT);']); + end + %% + SMALL.solver.time = cputime - start; + fprintf('Solver %s finished task in %2f seconds. \n', SMALL.solver.name, SMALL.solver.time); + SMALL.solver.solution = full(y); + SMALL.solver.reconstructed = SMALL.Problem.reconstruct(SMALL.solver.solution); +end + \ No newline at end of file