# HG changeset patch
# User idamnjanovic
# Date 1257438961 0
# Node ID 7750624e0c73c856ab7cfbae6160711e965135b7
# Parent  5181bee80bc17f1ebe496f114247b78e98a7f755


diff -r 5181bee80bc1 -r 7750624e0c73 README.txt
--- /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.
diff -r 5181bee80bc1 -r 7750624e0c73 SMALLboxSetup.m
--- /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
+

+  
diff -r 5181bee80bc1 -r 7750624e0c73 examples/SMALL_solver_test.m
--- /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
diff -r 5181bee80bc1 -r 7750624e0c73 examples/SMALL_solver_test_Audio.m
--- /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
diff -r 5181bee80bc1 -r 7750624e0c73 examples/private/add_dc.m
--- /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
+
+
+
diff -r 5181bee80bc1 -r 7750624e0c73 examples/private/addtocols.c
--- /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;
+}
diff -r 5181bee80bc1 -r 7750624e0c73 examples/private/addtocols.m
--- /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
diff -r 5181bee80bc1 -r 7750624e0c73 examples/private/addtocols.mexw32
Binary file examples/private/addtocols.mexw32 has changed
diff -r 5181bee80bc1 -r 7750624e0c73 examples/private/collincomb.c
--- /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;
+}
diff -r 5181bee80bc1 -r 7750624e0c73 examples/private/collincomb.m
--- /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
diff -r 5181bee80bc1 -r 7750624e0c73 examples/private/collincomb.mexw32
Binary file examples/private/collincomb.mexw32 has changed
diff -r 5181bee80bc1 -r 7750624e0c73 examples/private/countcover.m
--- /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
diff -r 5181bee80bc1 -r 7750624e0c73 examples/private/diag_ids.m
--- /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
diff -r 5181bee80bc1 -r 7750624e0c73 examples/private/dictdist.m
--- /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;
diff -r 5181bee80bc1 -r 7750624e0c73 examples/private/imnormalize.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);
diff -r 5181bee80bc1 -r 7750624e0c73 examples/private/iswhole.m
--- /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;
diff -r 5181bee80bc1 -r 7750624e0c73 examples/private/make.m
--- /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
diff -r 5181bee80bc1 -r 7750624e0c73 examples/private/normcols.m
--- /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)));
diff -r 5181bee80bc1 -r 7750624e0c73 examples/private/printf.m
--- /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
diff -r 5181bee80bc1 -r 7750624e0c73 examples/private/reggrid.m
--- /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
+
diff -r 5181bee80bc1 -r 7750624e0c73 examples/private/remove_dc.m
--- /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
+
diff -r 5181bee80bc1 -r 7750624e0c73 examples/private/rowlincomb.c
--- /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;
+}
diff -r 5181bee80bc1 -r 7750624e0c73 examples/private/rowlincomb.m
--- /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
diff -r 5181bee80bc1 -r 7750624e0c73 examples/private/rowlincomb.mexw32
Binary file examples/private/rowlincomb.mexw32 has changed
diff -r 5181bee80bc1 -r 7750624e0c73 examples/private/sampgrid.m
--- /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
+
diff -r 5181bee80bc1 -r 7750624e0c73 examples/private/secs2hms.m
--- /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);
diff -r 5181bee80bc1 -r 7750624e0c73 examples/private/spdiag.m
--- /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
diff -r 5181bee80bc1 -r 7750624e0c73 examples/private/timerclear.m
--- /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 = [];
diff -r 5181bee80bc1 -r 7750624e0c73 examples/private/timereta.m
--- /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
+
diff -r 5181bee80bc1 -r 7750624e0c73 examples/private/timerinit.m
--- /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
diff -r 5181bee80bc1 -r 7750624e0c73 solvers/SMALL_MP.m
--- /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;
+
+
diff -r 5181bee80bc1 -r 7750624e0c73 solvers/SMALL_cgp.m
--- /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;
+
+
+
diff -r 5181bee80bc1 -r 7750624e0c73 solvers/SMALL_chol.m
--- /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
diff -r 5181bee80bc1 -r 7750624e0c73 util/SL_A.m
--- /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
diff -r 5181bee80bc1 -r 7750624e0c73 util/SMALL_playAudio.m
--- /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
diff -r 5181bee80bc1 -r 7750624e0c73 util/SMALL_plot.m
--- /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
diff -r 5181bee80bc1 -r 7750624e0c73 util/SMALL_solve.m
--- /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