Mercurial > hg > smallbox
changeset 140:31d2864dfdd4 ivand_dev
Audio Impainting additional constraints with cvx added
| author | Ivan Damnjanovic lnx <ivan.damnjanovic@eecs.qmul.ac.uk> |
|---|---|
| date | Mon, 25 Jul 2011 17:27:05 +0100 |
| parents | 4bd6856a7128 |
| children | b8db4285e0ca |
| files | Problems/AudioDeclipping_reconstruct.m Problems/generateAudioDeclippingProblem.m SMALLboxSetup.m examples/AudioInpainting/Audio_Declipping_Example.m solvers/CVX_add_const_Audio_declipping.m solvers/SMALL_ompGabor/make.m solvers/SMALL_ompGabor/mexutils.c solvers/SMALL_ompGabor/mexutils.h solvers/SMALL_ompGabor/myblas.c solvers/SMALL_ompGabor/myblas.h solvers/SMALL_ompGabor/omp2Gabor.m solvers/SMALL_ompGabor/omp2mex.c solvers/SMALL_ompGabor/omp2mex.m solvers/SMALL_ompGabor/omp2mexGabor.c solvers/SMALL_ompGabor/omp2mexGabor.m solvers/SMALL_ompGabor/ompGabor.m solvers/SMALL_ompGabor/ompcore.c solvers/SMALL_ompGabor/ompcore.h solvers/SMALL_ompGabor/ompcoreGabor.c solvers/SMALL_ompGabor/ompcoreGabor.h solvers/SMALL_ompGabor/ompmex.c solvers/SMALL_ompGabor/ompmex.m solvers/SMALL_ompGabor/ompmexGabor.c solvers/SMALL_ompGabor/ompmexGabor.m solvers/SMALL_ompGabor/ompprof.c solvers/SMALL_ompGabor/ompprof.h solvers/SMALL_ompGabor/omputils.c solvers/SMALL_ompGabor/omputils.h solvers/SMALL_ompGabor/printf.m toolboxes/AudioInpaintingToolbox/Experiments/DeclippingExperiment/declipOneSoundExperiment.m util/SMALL_init_solver.m util/SMALL_solve.m util/ksvd utils/ompbox utils/make.m util/ksvd utils/ompbox utils/mexutils.c util/ksvd utils/ompbox utils/mexutils.h util/ksvd utils/ompbox utils/myblas.c util/ksvd utils/ompbox utils/myblas.h util/ksvd utils/ompbox utils/omp2Gabor.m util/ksvd utils/ompbox utils/omp2mex.c util/ksvd utils/ompbox utils/omp2mex.m util/ksvd utils/ompbox utils/omp2mexGabor.c util/ksvd utils/ompbox utils/omp2mexGabor.m util/ksvd utils/ompbox utils/ompGabor.m util/ksvd utils/ompbox utils/ompcore.c util/ksvd utils/ompbox utils/ompcore.h util/ksvd utils/ompbox utils/ompcoreGabor.c util/ksvd utils/ompbox utils/ompcoreGabor.h util/ksvd utils/ompbox utils/ompmex.c util/ksvd utils/ompbox utils/ompmex.m util/ksvd utils/ompbox utils/ompmexGabor.c util/ksvd utils/ompbox utils/ompmexGabor.m util/ksvd utils/ompbox utils/ompprof.c util/ksvd utils/ompbox utils/ompprof.h util/ksvd utils/ompbox utils/omputils.c util/ksvd utils/ompbox utils/omputils.h util/ksvd utils/ompbox utils/printf.m |
| diffstat | 56 files changed, 4187 insertions(+), 4021 deletions(-) [+] |
line wrap: on
line diff
--- a/Problems/AudioDeclipping_reconstruct.m Thu Jul 21 16:37:14 2011 +0100 +++ b/Problems/AudioDeclipping_reconstruct.m Mon Jul 25 17:27:05 2011 +0100 @@ -2,8 +2,15 @@ %% Audio declipping Problem reconstruction function % % This reconstruction function is using sparse representation y -% in dictionary Problem.A to reconstruct the patches of the denoised -% image. +% in dictionary Problem.A to reconstruct declipped audio. +% +% [1] I. Damnjanovic, M. E. P. Davies, and M. P. Plumbley "SMALLbox - an +% evaluation framework for sparse representations and dictionary +% learning algorithms," V. Vigneron et al. (Eds.): LVA/ICA 2010, +% Springer-Verlag, Berlin, Germany, LNCS 6365, pp. 418-425 +% [2] A. Adler, V. Emiya, M. G. Jafari, M. Elad, R. Gribonval, and M. D. +% Plumbley, “Audio Inpainting,” submitted to IEEE Trans. Audio, Speech, +% and Lang. Proc., 2011, http://hal.inria.fr/inria-00577079/en/. % % Centre for Digital Music, Queen Mary, University of London.
--- a/Problems/generateAudioDeclippingProblem.m Thu Jul 21 16:37:14 2011 +0100 +++ b/Problems/generateAudioDeclippingProblem.m Mon Jul 25 17:27:05 2011 +0100 @@ -1,26 +1,18 @@ function data = generateAudioDeclippingProblem(soundfile, clippingLevel, windowSize, overlap, wa, ws, wd, Dict_fun, redundancyFactor) %% Generate Audio Declipping Problem % -% CHANGE!!!generateAMT_Learning_Problem is a part of the SMALLbox and generates -% a problem that can be used for comparison of Dictionary Learning/Sparse -% Representation techniques in automatic music transcription scenario. -% The function prompts a user for an audio file (mid, wav, mat) reads it -% and generates a spectrogram given fft size (default nfft=4096), analysis -% window size (windowSize=2822), and analysis window overlap (overlap = -% 0.5). -% -% The output of the function is stucture with following fields: -% b - matrix with magnitudes of the spectrogram -% f - vector of frequencies at wihch spectrogram is computed -% windowSize - analysis window size -% overlap - analysis window overlap -% fs - sampling frequency -% m - number of frequenciy points in spectrogram -% n - number of time points in the spectrogram -% p - number of dictionary elements to be learned (eg 88 for piano) -% notesOriginal - notes of the original audio to be used for -% comparison (if midi of the original exists) -% name - name of the audio file to transcribe +% generateAudioDeclippingProblem is part of the SMALLbox [1] and generates +% Audio declipping is a problem proposed in Audio Inpaining Toolbox and +% in [2]. +% +% [1] I. Damnjanovic, M. E. P. Davies, and M. P. Plumbley "SMALLbox - an +% evaluation framework for sparse representations and dictionary +% learning algorithms," V. Vigneron et al. (Eds.): LVA/ICA 2010, +% Springer-Verlag, Berlin, Germany, LNCS 6365, pp. 418-425 +% [2] A. Adler, V. Emiya, M. G. Jafari, M. Elad, R. Gribonval, and M. D. +% Plumbley, “Audio Inpainting,” submitted to IEEE Trans. Audio, Speech, +% and Lang. Proc., 2011, http://hal.inria.fr/inria-00577079/en/. + % Centre for Digital Music, Queen Mary, University of London. % This file copyright 2011 Ivan Damnjanovic.
--- a/SMALLboxSetup.m Thu Jul 21 16:37:14 2011 +0100 +++ b/SMALLboxSetup.m Mon Jul 25 17:27:05 2011 +0100 @@ -274,7 +274,15 @@ end %% +%% KSVD utils setup +if ~(exist('addtocols')==3) + cd([SMALL_path,FS,'util',FS,'ksvd utils']); + make + cd(SMALL_path); +end + +%% if ~exist('ksvdver.m','file') fprintf('\n ******************************************************************'); fprintf('\n\n Initialising OMPbox and KSVDBox Setup'); @@ -300,9 +308,9 @@ 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); + 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'); @@ -392,13 +400,7 @@ fprintf('\n\n matlab_midi (http://www.kenschutte.com/midi/) is already installed'); end -%% KSVD utils setup -if ~(exist('addtocols')==3) - cd([SMALL_path,FS,'util',FS,'ksvd utils']); - make - cd(SMALL_path); -end %% RWT setup if ~(exist('mdwt')==3)
--- a/examples/AudioInpainting/Audio_Declipping_Example.m Thu Jul 21 16:37:14 2011 +0100 +++ b/examples/AudioInpainting/Audio_Declipping_Example.m Mon Jul 25 17:27:05 2011 +0100 @@ -1,11 +1,16 @@ %% Audio Declipping Example % -% CHANGE!!! This example is based on the experiment suggested by Professor Pierre -% Vandergheynst on the SMALL meeting in Villars. -% The idea behind is to use patches from source image as a dictionary in -% which we represent target image using matching pursuit algorithm. -% Calling Pierre_Problem function to get src image to be used as dictionary -% and target image to be represented using MP with 3 patches from source image +% Audio declipping is a problem proposed in Audio Inpaining Toolbox and +% in [2]. This is an example of solving the problem with fast omp using +% Gabor dictionary and ompGabor implemented in SMALLbox [1]. +% +% [1] I. Damnjanovic, M. E. P. Davies, and M. P. Plumbley "SMALLbox - an +% evaluation framework for sparse representations and dictionary +% learning algorithms," V. Vigneron et al. (Eds.): LVA/ICA 2010, +% Springer-Verlag, Berlin, Germany, LNCS 6365, pp. 418-425 +% [2] A. Adler, V. Emiya, M. G. Jafari, M. Elad, R. Gribonval, and M. D. +% Plumbley, “Audio Inpainting,” submitted to IEEE Trans. Audio, Speech, +% and Lang. Proc., 2011, http://hal.inria.fr/inria-00577079/en/. % % Centre for Digital Music, Queen Mary, University of London. @@ -20,88 +25,123 @@ %% clear all; +% Defining the solvers to test in Audio declipping scenario + +% First solver omp2 - DCT+DST dictionary with no additional constraints + +SMALL.solver(1) = SMALL_init_solver('ompbox', 'omp2', '', 0); +SMALL.solver(1).add_constraints = 0; + +% Second solver omp2 - DCT+DST dictionary with additional constraints + +SMALL.solver(2) = SMALL_init_solver('ompbox', 'omp2', '', 0); +SMALL.solver(2).add_constraints = 1; + +% Third solver omp2 - Gabor dictionary with no additional constraints + +SMALL.solver(3) = SMALL_init_solver('ompbox', 'omp2Gabor', '', 0); +SMALL.solver(3).add_constraints = 0; + +% Fourth solver omp2- Gabor dictionary with no additional constraints + +SMALL.solver(4) = SMALL_init_solver('ompbox', 'omp2Gabor', '', 0); +SMALL.solver(4).add_constraints = 1; % Defining the Problem structure SMALL.Problem = generateAudioDeclippingProblem('male01_8kHz.wav', 0.6, 256, 0.5, @wRect, @wSine, @wRect, @Gabor_Dictionary, 2); -% % Show original image and image that is used as a dictionary -% figure('Name', 'Original and Dictionary Image'); -% -% subplot(1,2,1); imagesc(SMALL.Problem.imageTrg/SMALL.Problem.maxval); -% title('Original Image');colormap(gray);axis off; axis image; -% subplot(1,2,2); imagesc(SMALL.Problem.imageSrc/SMALL.Problem.maxval); -% title('Dictionary image:');colormap(gray);axis off; axis image; -time=0; -error2=0.001^2; -coeffFrames = zeros(SMALL.Problem.p, SMALL.Problem.n); - -for i=1:SMALL.Problem.n +for idxSolver = 1:4 - idx = find(SMALL.Problem.M(:,i)); - if size(idx,1)==SMALL.Problem.m - continue + fprintf('\nStarting Audio Declipping of %s... \n', SMALL.Problem.name); + fprintf('\nClipping level %s... \n', SMALL.Problem.clippingLevel); + + start=cputime; + tStart=tic; + + error2=0.001^2; + coeffFrames = zeros(SMALL.Problem.p, SMALL.Problem.n); + + + + for i = 1:SMALL.Problem.n + + idx = find(SMALL.Problem.M(:,i)); + if size(idx,1)==SMALL.Problem.m + continue + end + Dict = SMALL.Problem.B(idx,:); + wDict = 1./sqrt(diag(Dict'*Dict)); + + SMALL.Problem.A = Dict*diag(wDict); + + SMALL.Problem.b1 = SMALL.Problem.b(idx,i); + + + + % set solver parameters + + SMALL.solver(idxSolver).param=struct(... + 'epsilon', error2*size(idx,1),... + 'maxatoms', 128, ... + 'profile', 'off'); + + % Find solution + + SMALL.solver(idxSolver)=SMALL_solve(SMALL.Problem, SMALL.solver(idxSolver)); + + % Refine solution with additional constraints for declipping scenario + + if (SMALL.solver(idxSolver).add_constraints) + SMALL.solver(idxSolver)=CVX_add_const_Audio_declipping(... + SMALL.Problem, SMALL.solver(idxSolver), i); + end + + %% + coeffFrames(:,i) = wDict .* SMALL.solver(idxSolver).solution; + + end - Dict = SMALL.Problem.B(idx,:); - wDict = 1./sqrt(diag(Dict'*Dict)); - SMALL.Problem.A = Dict*diag(wDict); + %% Set reconstruction function - SMALL.Problem.b1 = SMALL.Problem.b(idx,i); + SMALL.Problem.reconstruct=@(x) AudioDeclipping_reconstruct(x, SMALL.Problem); + reconstructed=SMALL.Problem.reconstruct(coeffFrames); + SMALL.Problem = rmfield(SMALL.Problem, 'reconstruct'); + tElapsed=toc(tStart); + SMALL.solver(idxSolver).time = cputime - start; + fprintf('Solver %s finished task in %2f seconds (cpu time). \n', ... + SMALL.solver(idxSolver).name, SMALL.solver(idxSolver).time); + fprintf('Solver %s finished task in %2f seconds (tic-toc time). \n', ... + SMALL.solver(idxSolver).name, tElapsed); - % Defining the parameters sparse representation - SMALL.solver=SMALL_init_solver; - SMALL.solver.toolbox='ompbox'; - SMALL.solver.name='omp2Gabor'; - SMALL.solver.param=struct(... - 'epsilon', error2*size(idx,1),... - 'maxatoms', 128); + %% Plot results - % Find solution + xClipped = SMALL.Problem.clipped; + xClean = SMALL.Problem.original; + xEst1 = reconstructed.audioAllSamples; + xEst2 = reconstructed.audioOnlyClipped; + fs = SMALL.Problem.fs; - SMALL.solver=SMALL_solve(SMALL.Problem, SMALL.solver); - - - coeffFrames(:,i) = wDict .* SMALL.solver.solution; - time = time + SMALL.solver.time; - - - + figure + hold on + plot(xClipped,'r') + plot(xClean) + plot(xEst2,'--g') + plot([1;length(xClipped)],[1;1]*[-1,1]*max(abs(xClipped)),':r') + legend('Clipped','True solution','Estimate') end -%% Set reconstruction function - -SMALL.Problem.reconstruct=@(x) AudioDeclipping_reconstruct(x, SMALL.Problem); -reconstructed=SMALL.Problem.reconstruct(coeffFrames); - - - -%% Plot results - -xClipped = SMALL.Problem.clipped; -xClean = SMALL.Problem.original; -xEst1 = reconstructed.audioAllSamples; -xEst2 = reconstructed.audioOnlyClipped; -fs = SMALL.Problem.fs; - -figure -hold on -plot(xClipped,'r') -plot(xClean) -plot(xEst2,'--g') -plot([1;length(xClipped)],[1;1]*[-1,1]*max(abs(xClipped)),':r') -legend('Clipped','True solution','Estimate') - -% Normalized and save sounds -normX = 1.1*max(abs([xEst1(:);xEst2(:);xClean(:)])); -L = min([length(xEst2),length(xEst1),length(xClean),length(xEst1),length(xClipped)]); -xEst1 = xEst1(1:L)/normX; -xEst2 = xEst2(1:L)/normX; -xClipped = xClipped(1:L)/normX; -xClean = xClean(1:L)/normX; +% % Normalized and save sounds +% normX = 1.1*max(abs([xEst1(:);xEst2(:);xClean(:)])); +% L = min([length(xEst2),length(xEst1),length(xClean),length(xEst1),length(xClipped)]); +% xEst1 = xEst1(1:L)/normX; +% xEst2 = xEst2(1:L)/normX; +% xClipped = xClipped(1:L)/normX; +% xClean = xClean(1:L)/normX; % wavwrite(xEst1,fs,[expParam.destDir 'xEst1.wav']); % wavwrite(xEst2,fs,[expParam.destDir 'xEst2.wav']); % wavwrite(xClipped,fs,[expParam.destDir 'xClipped.wav']);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/solvers/CVX_add_const_Audio_declipping.m Mon Jul 25 17:27:05 2011 +0100 @@ -0,0 +1,85 @@ +function solver=CVX_add_const_Audio_declipping(Problem, solver, idxFrame) + +%% CVX additional constrains + % Clipping level: take the analysis window into account + % clipping level detection + + signal = Problem.b1; + signalFull = Problem.b(:,idxFrame); + Dict = Problem.A; + DictFull = Problem.B; + Clipped = ~Problem.M(:,idxFrame); + W=1./sqrt(diag(Dict'*Dict)); + + idxCoeff = find(solver.solution~=0); + + solution = solver.solution; + + wa = Problem.wa(Problem.windowSize); % analysis window + + + clippingLevelEst = max(abs(signal./wa(~Clipped)')); + + idxPos = find(signalFull>=0 & Clipped); + idxNeg = find(signalFull<0 & Clipped); + + DictPos=DictFull(idxPos,:); + DictNeg=DictFull(idxNeg,:); + + + wa_pos = wa(idxPos); + wa_neg = wa(idxNeg); + b_ineq_pos = wa_pos(:)*clippingLevelEst; + b_ineq_neg = -wa_neg(:)*clippingLevelEst; + if isfield(Problem,'Upper_Limit') && ~isempty(Problem.Upper_Limit) + b_ineq_pos_upper_limit = wa_pos(:)*Problem.Upper_Limit*clippingLevelEst; + b_ineq_neg_upper_limit = -wa_neg(:)*Problem.Upper_Limit*clippingLevelEst; + else + b_ineq_pos_upper_limit = Inf; + b_ineq_neg_upper_limit = -Inf; + end + + + j = length(idxCoeff); + + if isinf(b_ineq_pos_upper_limit) + %% CVX code + cvx_begin + cvx_quiet(true) + variable solution(j) + minimize(norm(Dict(:,idxCoeff)*solution-signal)) + subject to + DictPos(:,idxCoeff)*(W(idxCoeff).*solution) >= b_ineq_pos + DictNeg(:,idxCoeff)*(W(idxCoeff).*solution) <= b_ineq_neg + cvx_end + if cvx_optval>1e3 + cvx_begin + cvx_quiet(true) + variable solution(j) + minimize(norm(Dict(:,idxCoeff)*solution-xObs)) + cvx_end + end + else + %% CVX code + cvx_begin + cvx_quiet(true) + variable solution(j) + minimize(norm(Dict(:,idxCoeff)*solution-signal)) + subject to + DictPos(:,idxCoeff)*(W(idxCoeff).*solution) >= b_ineq_pos + DictNeg(:,idxCoeff)*(W(idxCoeff).*solution) <= b_ineq_neg + DictPos(:,idxCoeff)*(W(idxCoeff).*solution) <= b_ineq_pos_upper_limit + DictNeg(:,idxCoeff)*(W(idxCoeff).*solution) >= b_ineq_neg_upper_limit + cvx_end + if cvx_optval>1e3 + cvx_begin + cvx_quiet(true) + variable solution(j) + minimize(norm(Dict(:,idxCoeff)*solution-xObs)) + cvx_end + end + end + + solver.solution(idxCoeff) = solution; + + \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/solvers/SMALL_ompGabor/make.m Mon Jul 25 17:27:05 2011 +0100 @@ -0,0 +1,41 @@ +function make +%MAKE Build the OMPBox package. +% MAKE compiles all OMPBox MEX functions, using Matlab's default MEX +% compiler. If the MEX compiler has not been set-up before, please run +% +% mex -setup +% +% before using this MAKE file. + +% Ron Rubinstein +% Computer Science Department +% Technion, Haifa 32000 Israel +% ronrubin@cs +% +% August 2009 + + +% detect platform + +compstr = computer; +is64bit = strcmp(compstr(end-1:end),'64'); + + +% compilation parameters + +compile_params = cell(0); +if (is64bit) + compile_params{1} = '-largeArrayDims'; +end + + +% Compile files % + +ompsources = {'mexutils.c','ompcoreGabor.c','omputils.c','myblas.c','ompprof.c'}; + +disp('Compiling ompmex...'); +mex('ompmexGabor.c', ompsources{:},compile_params{:}); + +disp('Compiling omp2mex...'); +mex('omp2mexGabor.c',ompsources{:},compile_params{:}); +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/solvers/SMALL_ompGabor/mexutils.c Mon Jul 25 17:27:05 2011 +0100 @@ -0,0 +1,79 @@ +/************************************************************************** + * + * File name: mexutils.c + * + * Ron Rubinstein + * Computer Science Department + * Technion, Haifa 32000 Israel + * ronrubin@cs + * + * Last Updated: 15.8.2009 + * + *************************************************************************/ + +#include "mexutils.h" +#include <math.h> + + + +/* verify that the mxArray contains a double matrix */ + +void checkmatrix(const mxArray *param, char *fname, char *pname) +{ + char errmsg[100]; + sprintf(errmsg, "%.15s requires that %.25s be a double matrix.", fname, pname); + if (!mxIsDouble(param) || mxIsComplex(param) || mxGetNumberOfDimensions(param)>2) { + mexErrMsgTxt(errmsg); + } +} + + +/* verify that the mxArray contains a 1-D double vector */ + +void checkvector(const mxArray *param, char *fname, char *pname) +{ + char errmsg[100]; + sprintf(errmsg, "%.15s requires that %.25s be a double vector.", fname, pname); + if (!mxIsDouble(param) || mxIsComplex(param) || mxGetNumberOfDimensions(param)>2 || (mxGetM(param)!=1 && mxGetN(param)!=1)) { + mexErrMsgTxt(errmsg); + } +} + + +/* verify that the mxArray contains a double scalar */ + +void checkscalar(const mxArray *param, char *fname, char *pname) +{ + char errmsg[100]; + sprintf(errmsg, "%.15s requires that %.25s be a double scalar.", fname, pname); + if (!mxIsDouble(param) || mxIsComplex(param) || mxGetNumberOfDimensions(param)>2 || + mxGetM(param)!=1 || mxGetN(param)!=1) + { + mexErrMsgTxt(errmsg); + } +} + + +/* verify that the mxArray contains a sparse matrix */ + +void checksparse(const mxArray *param, char *fname, char *pname) +{ + char errmsg[100]; + sprintf(errmsg, "%.15s requires that %.25s be sparse.", fname, pname); + if (!mxIsSparse(param)) { + mexErrMsgTxt(errmsg); + } +} + + +/* verify that the mxArray contains a 1-dimensional cell array */ + +void checkcell_1d(const mxArray *param, char *fname, char *pname) +{ + char errmsg[100]; + sprintf(errmsg, "%.15s requires that %.25s be a 1-D cell array.", fname, pname); + if (!mxIsCell(param) || mxGetNumberOfDimensions(param)>2 || (mxGetM(param)!=1 && mxGetN(param)!=1)) { + mexErrMsgTxt(errmsg); + } +} +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/solvers/SMALL_ompGabor/mexutils.h Mon Jul 25 17:27:05 2011 +0100 @@ -0,0 +1,103 @@ +/************************************************************************** + * + * File name: mexutils.h + * + * Ron Rubinstein + * Computer Science Department + * Technion, Haifa 32000 Israel + * ronrubin@cs + * + * Last Updated: 18.8.2009 + * + * Utility functions for MEX files. + * + *************************************************************************/ + + +#ifndef __MEX_UTILS_H__ +#define __MEX_UTILS_H__ + +#include "mex.h" + + + +/************************************************************************** + * Function checkmatrix: + * + * Verify that the specified mxArray is real, of type double, and has + * no more than two dimensions. If not, an error message is printed + * and the mex file terminates. + * + * Parameters: + * param - the mxArray to be checked + * fname - the name of the function where the error occured (15 characters or less) + * pname - the name of the parameter (25 characters or less) + * + **************************************************************************/ +void checkmatrix(const mxArray *param, char *fname, char *pname); + + +/************************************************************************** + * Function checkvector: + * + * Verify that the specified mxArray is 1-D, real, and of type double. The + * vector may be a column or row vector. Otherwise, an error message is + * printed and the mex file terminates. + * + * Parameters: + * param - the mxArray to be checked + * fname - the name of the function where the error occured (15 characters or less) + * pname - the name of the parameter (25 characters or less) + * + **************************************************************************/ +void checkvector(const mxArray *param, char *fname, char *pname); + + +/************************************************************************** + * Function checkscalar: + * + * Verify that the specified mxArray represents a real double scalar value. + * If not, an error message is printed and the mex file terminates. + * + * Parameters: + * param - the mxArray to be checked + * fname - the name of the function where the error occured (15 characters or less) + * pname - the name of the parameter (25 characters or less) + * + **************************************************************************/ +void checkscalar(const mxArray *param, char *fname, char *pname); + + +/************************************************************************** + * Function checksparse: + * + * Verify that the specified mxArray contains a sparse matrix. If not, + * an error message is printed and the mex file terminates. + * + * Parameters: + * param - the mxArray to be checked + * fname - the name of the function where the error occured (15 characters or less) + * pname - the name of the parameter (25 characters or less) + * + **************************************************************************/ +void checksparse(const mxArray *param, char *fname, char *pname); + + +/************************************************************************** + * Function checkcell_1d: + * + * Verify that the specified mxArray is a 1-D cell array. The cell array + * may be arranged as either a column or a row. If not, an error message + * is printed and the mex file terminates. + * + * Parameters: + * param - the mxArray to be checked + * fname - the name of the function where the error occured (15 characters or less) + * pname - the name of the parameter (25 characters or less) + * + **************************************************************************/ +void checkcell_1d(const mxArray *param, char *fname, char *pname); + + +#endif +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/solvers/SMALL_ompGabor/myblas.c Mon Jul 25 17:27:05 2011 +0100 @@ -0,0 +1,673 @@ +/************************************************************************** + * + * File name: myblas.c + * + * Ron Rubinstein + * Computer Science Department + * Technion, Haifa 32000 Israel + * ronrubin@cs + * + * Version: 1.1 + * Last updated: 13.8.2009 + * + *************************************************************************/ + + +#include "myblas.h" +#include <ctype.h> + + +/* find maximum of absolute values */ + +mwIndex maxabs(double c[], mwSize m) +{ + mwIndex maxid=0, k; + double absval, maxval = SQR(*c); /* use square which is quicker than absolute value */ + + for (k=1; k<m; ++k) { + absval = SQR(c[k]); + if (absval > maxval) { + maxval = absval; + maxid = k; + } + } + return maxid; +} + + +/* compute y := alpha*x + y */ + +void vec_sum(double alpha, double x[], double y[], mwSize n) +{ + mwIndex i; + + for (i=0; i<n; ++i) { + y[i] += alpha*x[i]; + } +} + +/* compute y := alpha*x .* y */ + +void vec_smult(double alpha, double x[], double y[], mwSize n) +{ + mwIndex i; + + for (i=0; i<n; ++i) { + y[i] *= alpha*x[i]; + } +} + +/* compute y := alpha*A*x */ + +void mat_vec(double alpha, double A[], double x[], double y[], mwSize n, mwSize m) +{ + mwIndex i, j, i_n; + double *Ax; + + Ax = mxCalloc(n,sizeof(double)); + + for (i=0; i<m; ++i) { + i_n = i*n; + for (j=0; j<n; ++j) { + Ax[j] += A[i_n+j] * x[i]; + } + } + + for (j=0; j<n; ++j) { + y[j] = alpha*Ax[j]; + } + + mxFree(Ax); +} + + +/* compute y := alpha*A'*x */ + +void matT_vec(double alpha, double A[], double x[], double y[], mwSize n, mwSize m) +{ + mwIndex i, j, n_i; + double sum0, sum1, sum2, sum3; + + for (j=0; j<m; ++j) { + y[j] = 0; + } + + /* use loop unrolling to accelerate computation */ + + for (i=0; i<m; ++i) { + n_i = n*i; + sum0 = sum1 = sum2 = sum3 = 0; + for (j=0; j+4<n; j+=4) { + sum0 += A[n_i+j]*x[j]; + sum1 += A[n_i+j+1]*x[j+1]; + sum2 += A[n_i+j+2]*x[j+2]; + sum3 += A[n_i+j+3]*x[j+3]; + } + y[i] += alpha * ((sum0 + sum1) + (sum2 + sum3)); + while (j<n) { + y[i] += alpha*A[n_i+j]*x[j]; + j++; + } + } +} + + +/* compute y := alpha*A*x */ + +void mat_sp_vec(double alpha, double pr[], mwIndex ir[], mwIndex jc[], double x[], double y[], mwSize n, mwSize m) +{ + + mwIndex i, j, j1, j2; + + for (i=0; i<n; ++i) { + y[i] = 0; + } + + j2 = jc[0]; + for (i=0; i<m; ++i) { + j1 = j2; j2 = jc[i+1]; + for (j=j1; j<j2; ++j) { + y[ir[j]] += alpha * pr[j] * x[i]; + } + } + +} + + +/* compute y := alpha*A'*x */ + +void matT_sp_vec(double alpha, double pr[], mwIndex ir[], mwIndex jc[], double x[], double y[], mwSize n, mwSize m) +{ + + mwIndex i, j, j1, j2; + + for (i=0; i<m; ++i) { + y[i] = 0; + } + + j2 = jc[0]; + for (i=0; i<m; ++i) { + j1 = j2; j2 = jc[i+1]; + for (j=j1; j<j2; ++j) { + y[i] += alpha * pr[j] * x[ir[j]]; + } + } + +} + + +/* compute y := alpha*A*x */ + +void mat_vec_sp(double alpha, double A[], double pr[], mwIndex ir[], mwIndex jc[], double y[], mwSize n, mwSize m) +{ + + mwIndex i, j, j_n, k, kend; + + for (i=0; i<n; ++i) { + y[i] = 0; + } + + kend = jc[1]; + if (kend==0) { /* x is empty */ + return; + } + + for (k=0; k<kend; ++k) { + j = ir[k]; + j_n = j*n; + for (i=0; i<n; ++i) { + y[i] += alpha * A[i+j_n] * pr[k]; + } + } + +} + + +/* compute y := alpha*A'*x */ + +void matT_vec_sp(double alpha, double A[], double pr[], mwIndex ir[], mwIndex jc[], double y[], mwSize n, mwSize m) +{ + + mwIndex i, j, j_n, k, kend; + + for (i=0; i<m; ++i) { + y[i] = 0; + } + + kend = jc[1]; + if (kend==0) { /* x is empty */ + return; + } + + for (j=0; j<m; ++j) { + j_n = j*n; + for (k=0; k<kend; ++k) { + i = ir[k]; + y[j] += alpha * A[i+j_n] * pr[k]; + } + } + +} + + +/* compute y := alpha*A*x */ + +void mat_sp_vec_sp(double alpha, double pr[], mwIndex ir[], mwIndex jc[], double prx[], mwIndex irx[], mwIndex jcx[], double y[], mwSize n, mwSize m) +{ + + mwIndex i, j, k, kend, j1, j2; + + for (i=0; i<n; ++i) { + y[i] = 0; + } + + kend = jcx[1]; + if (kend==0) { /* x is empty */ + return; + } + + for (k=0; k<kend; ++k) { + i = irx[k]; + j1 = jc[i]; j2 = jc[i+1]; + for (j=j1; j<j2; ++j) { + y[ir[j]] += alpha * pr[j] * prx[k]; + } + } + +} + + +/* compute y := alpha*A'*x */ + +void matT_sp_vec_sp(double alpha, double pr[], mwIndex ir[], mwIndex jc[], double prx[], mwIndex irx[], mwIndex jcx[], double y[], mwSize n, mwSize m) +{ + + mwIndex i, j, k, jend, kend, jadd, kadd, delta; + + for (i=0; i<m; ++i) { + y[i] = 0; + } + + kend = jcx[1]; + if (kend==0) { /* x is empty */ + return; + } + + for (i=0; i<m; ++i) { + j = jc[i]; + jend = jc[i+1]; + k = 0; + while (j<jend && k<kend) { + + delta = ir[j] - irx[k]; + + if (delta) { /* if indices differ - increment the smaller one */ + jadd = delta<0; + kadd = 1-jadd; + j += jadd; + k += kadd; + } + + else { /* indices are equal - add to result and increment both */ + y[i] += alpha * pr[j] * prx[k]; + j++; k++; + } + } + } + +} + + +/* matrix-matrix multiplication */ + +void mat_mat(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k) +{ + mwIndex i1, i2, i3, iX, iA, i2_n; + double b; + + for (i1=0; i1<n*k; i1++) { + X[i1] = 0; + } + + for (i2=0; i2<m; ++i2) { + i2_n = i2*n; + iX = 0; + for (i3=0; i3<k; ++i3) { + iA = i2_n; + b = B[i2+i3*m]; + for (i1=0; i1<n; ++i1) { + X[iX++] += A[iA++]*b; + } + } + } + + for (i1=0; i1<n*k; i1++) { + X[i1] *= alpha; + } +} + + +/* matrix-transpose-matrix multiplication */ + +void matT_mat(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k) +{ + mwIndex i1, i2, i3, iX, iA, i2_n; + double *x, sum0, sum1, sum2, sum3; + + for (i2=0; i2<m; ++i2) { + for (i3=0; i3<k; ++i3) { + sum0 = sum1 = sum2 = sum3 = 0; + for (i1=0; i1+4<n; i1+=4) { + sum0 += A[i1+0+i2*n]*B[i1+0+i3*n]; + sum1 += A[i1+1+i2*n]*B[i1+1+i3*n]; + sum2 += A[i1+2+i2*n]*B[i1+2+i3*n]; + sum3 += A[i1+3+i2*n]*B[i1+3+i3*n]; + } + X[i2+i3*m] = (sum0+sum1) + (sum2+sum3); + while(i1<n) { + X[i2+i3*m] += A[i1+i2*n]*B[i1+i3*n]; + i1++; + } + } + } + + for (i1=0; i1<m*k; i1++) { + X[i1] *= alpha; + } +} + + +/* tensor-matrix product */ + +void tens_mat(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k, mwSize l) +{ + mwIndex i1, i2, i3, i4, i2_n, nml; + double b; + + nml = n*m*l; + for (i1=0; i1<nml; ++i1) { + X[i1] = 0; + } + + for (i2=0; i2<m; ++i2) { + i2_n = i2*n; + for (i3=0; i3<k; ++i3) { + for (i4=0; i4<l; ++i4) { + b = B[i4+i3*l]; + for (i1=0; i1<n; ++i1) { + X[i1 + i2_n + i4*n*m] += A[i1 + i2_n + i3*n*m] * b; + } + } + } + } + + for (i1=0; i1<nml; ++i1) { + X[i1] *= alpha; + } +} + + +/* tensor-matrix-transpose product */ + +void tens_matT(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k, mwSize l) +{ + mwIndex i1, i2, i3, i4, i2_n, nml; + double b; + + nml = n*m*l; + for (i1=0; i1<nml; ++i1) { + X[i1] = 0; + } + + for (i2=0; i2<m; ++i2) { + i2_n = i2*n; + for (i4=0; i4<l; ++i4) { + for (i3=0; i3<k; ++i3) { + b = B[i3+i4*k]; + for (i1=0; i1<n; ++i1) { + X[i1 + i2_n + i4*n*m] += A[i1 + i2_n + i3*n*m] * b; + } + } + } + } + + for (i1=0; i1<nml; ++i1) { + X[i1] *= alpha; + } +} + + +/* dot product */ + +double dotprod(double a[], double b[], mwSize n) +{ + double sum = 0; + mwIndex i; + for (i=0; i<n; ++i) + sum += a[i]*b[i]; + return sum; +} + + +/* find maximum of vector */ + +mwIndex maxpos(double c[], mwSize m) +{ + mwIndex maxid=0, k; + double val, maxval = *c; + + for (k=1; k<m; ++k) { + val = c[k]; + if (val > maxval) { + maxval = val; + maxid = k; + } + } + return maxid; +} + + +/* solve L*x = b */ + +void backsubst_L(double L[], double b[], double x[], mwSize n, mwSize k) +{ + mwIndex i, j; + double rhs; + + for (i=0; i<k; ++i) { + rhs = b[i]; + for (j=0; j<i; ++j) { + rhs -= L[j*n+i]*x[j]; + } + x[i] = rhs/L[i*n+i]; + } +} + + +/* solve L'*x = b */ + +void backsubst_Lt(double L[], double b[], double x[], mwSize n, mwSize k) +{ + mwIndex i, j; + double rhs; + + for (i=k; i>=1; --i) { + rhs = b[i-1]; + for (j=i; j<k; ++j) { + rhs -= L[(i-1)*n+j]*x[j]; + } + x[i-1] = rhs/L[(i-1)*n+i-1]; + } +} + + +/* solve U*x = b */ + +void backsubst_U(double U[], double b[], double x[], mwSize n, mwSize k) +{ + mwIndex i, j; + double rhs; + + for (i=k; i>=1; --i) { + rhs = b[i-1]; + for (j=i; j<k; ++j) { + rhs -= U[j*n+i-1]*x[j]; + } + x[i-1] = rhs/U[(i-1)*n+i-1]; + } +} + + +/* solve U'*x = b */ + +void backsubst_Ut(double U[], double b[], double x[], mwSize n, mwSize k) +{ + mwIndex i, j; + double rhs; + + for (i=0; i<k; ++i) { + rhs = b[i]; + for (j=0; j<i; ++j) { + rhs -= U[i*n+j]*x[j]; + } + x[i] = rhs/U[i*n+i]; + } +} + + +/* back substitution solver */ + +void backsubst(char ul, double A[], double b[], double x[], mwSize n, mwSize k) +{ + if (tolower(ul) == 'u') { + backsubst_U(A, b, x, n, k); + } + else if (tolower(ul) == 'l') { + backsubst_L(A, b, x, n, k); + } + else { + mexErrMsgTxt("Invalid triangular matrix type: must be ''U'' or ''L''"); + } +} + + +/* solve equation set using cholesky decomposition */ + +void cholsolve(char ul, double A[], double b[], double x[], mwSize n, mwSize k) +{ + double *tmp; + + tmp = mxMalloc(k*sizeof(double)); + + if (tolower(ul) == 'l') { + backsubst_L(A, b, tmp, n, k); + backsubst_Lt(A, tmp, x, n, k); + } + else if (tolower(ul) == 'u') { + backsubst_Ut(A, b, tmp, n, k); + backsubst_U(A, tmp, x, n, k); + } + else { + mexErrMsgTxt("Invalid triangular matrix type: must be either ''U'' or ''L''"); + } + + mxFree(tmp); +} + + +/* perform a permutation assignment y := x(ind(1:k)) */ + +void vec_assign(double y[], double x[], mwIndex ind[], mwSize k) +{ + mwIndex i; + + for (i=0; i<k; ++i) + y[i] = x[ind[i]]; +} + + +/* matrix transpose */ + +void transpose(double X[], double Y[], mwSize n, mwSize m) +{ + mwIndex i, j, i_m, j_n; + + if (n<m) { + for (j=0; j<m; ++j) { + j_n = j*n; + for (i=0; i<n; ++i) { + Y[j+i*m] = X[i+j_n]; + } + } + } + else { + for (i=0; i<n; ++i) { + i_m = i*m; + for (j=0; j<m; ++j) { + Y[j+i_m] = X[i+j*n]; + } + } + } +} + + +/* print contents of matrix */ + +void printmat(double A[], int n, int m, char* matname) +{ + int i, j; + mexPrintf("\n%s = \n\n", matname); + + if (n*m==0) { + mexPrintf(" Empty matrix: %d-by-%d\n\n", n, m); + return; + } + + for (i=0; i<n; ++i) { + for (j=0; j<m; ++j) + mexPrintf(" %lf", A[j*n+i]); + mexPrintf("\n"); + } + mexPrintf("\n"); +} + + +/* print contents of sparse matrix */ + +void printspmat(mxArray *a, char* matname) +{ + mwIndex *aJc = mxGetJc(a); + mwIndex *aIr = mxGetIr(a); + double *aPr = mxGetPr(a); + + int i; + + mexPrintf("\n%s = \n\n", matname); + + for (i=0; i<aJc[1]; ++i) + printf(" (%d,1) = %lf\n", aIr[i]+1,aPr[i]); + + mexPrintf("\n"); +} + + + +/* matrix multiplication using Winograd's algorithm */ + +/* +void mat_mat2(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k) +{ + + mwIndex i1, i2, i3, iX, iA, i2_n; + double b, *AA, *BB; + + AA = mxCalloc(n,sizeof(double)); + BB = mxCalloc(k,sizeof(double)); + + for (i1=0; i1<n*k; i1++) { + X[i1] = 0; + } + + for (i1=0; i1<n; ++i1) { + for (i2=0; i2<m/2; ++i2) { + AA[i1] += A[i1+2*i2*n]*A[i1+(2*i2+1)*n]; + } + } + + for (i2=0; i2<k; ++i2) { + for (i1=0; i1<m/2; ++i1) { + BB[i2] += B[2*i1+i2*m]*B[2*i1+1+i2*m]; + } + } + + for (i2=0; i2<k; ++i2) { + for (i3=0; i3<m/2; ++i3) { + for (i1=0; i1<n; ++i1) { + X[i1+i2*n] += (A[i1+(2*i3)*n]+B[2*i3+1+i2*m])*(A[i1+(2*i3+1)*n]+B[2*i3+i2*m]); + } + } + } + + if (m%2) { + for (i2=0; i2<k; ++i2) { + for (i1=0; i1<n; ++i1) { + X[i1+i2*n] += A[i1+(m-1)*n]*B[m-1+i2*m]; + } + } + } + + for (i2=0; i2<k; ++i2) { + for (i1=0; i1<n; ++i1) { + X[i1+i2*n] -= (AA[i1] + BB[i2]); + X[i1+i2*n] *= alpha; + } + } + + mxFree(AA); + mxFree(BB); +} +*/ + + + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/solvers/SMALL_ompGabor/myblas.h Mon Jul 25 17:27:05 2011 +0100 @@ -0,0 +1,511 @@ +/************************************************************************** + * + * File name: myblas.h + * + * Ron Rubinstein + * Computer Science Department + * Technion, Haifa 32000 Israel + * ronrubin@cs + * + * Version: 1.1 + * Last updated: 17.8.2009 + * + * A collection of basic linear algebra functions, in the spirit of the + * BLAS/LAPACK libraries. + * + *************************************************************************/ + + + +#ifndef __MY_BLAS_H__ +#define __MY_BLAS_H__ + + +#include "mex.h" +#include <math.h> + + + +/************************************************************************** + * Squared value. + **************************************************************************/ +#define SQR(X) ((X)*(X)) + + + +/************************************************************************** + * Matrix-vector multiplication. + * + * Computes an operation of the form: + * + * y := alpha*A*x + * + * Parameters: + * A - matrix of size n X m + * x - vector of length m + * y - output vector of length n + * alpha - real constant + * n, m - dimensions of A + * + * Note: This function re-writes the contents of y. + * + **************************************************************************/ +void mat_vec(double alpha, double A[], double x[], double y[], mwSize n, mwSize m); + + + +/************************************************************************** + * Matrix-transpose-vector multiplication. + * + * Computes an operation of the form: + * + * y := alpha*A'*x + * + * Parameters: + * A - matrix of size n X m + * x - vector of length n + * y - output vector of length m + * alpha - real constant + * n, m - dimensions of A + * + * Note: This function re-writes the contents of y. + * + **************************************************************************/ +void matT_vec(double alpha, double A[], double x[], double y[], mwSize n, mwSize m); + + + +/************************************************************************** + * Sparse-matrix-vector multiplication. + * + * Computes an operation of the form: + * + * y := alpha*A*x + * + * where A is a sparse matrix. + * + * Parameters: + * pr,ir,jc - sparse representation of the matrix A, of size n x m + * x - vector of length m + * y - output vector of length n + * alpha - real constant + * n, m - dimensions of A + * + * Note: This function re-writes the contents of y. + * + **************************************************************************/ +void mat_sp_vec(double alpha, double pr[], mwIndex ir[], mwIndex jc[], double x[], double y[], mwSize n, mwSize m); + + + +/************************************************************************** + * Sparse-matrix-transpose-vector multiplication. + * + * Computes an operation of the form: + * + * y := alpha*A'*x + * + * where A is a sparse matrix. + * + * Parameters: + * pr,ir,jc - sparse representation of the matrix A, of size n x m + * x - vector of length m + * y - output vector of length n + * alpha - real constant + * n, m - dimensions of A + * + * Note: This function re-writes the contents of y. + * + **************************************************************************/ +void matT_sp_vec(double alpha, double pr[], mwIndex ir[], mwIndex jc[], double x[], double y[], mwSize n, mwSize m); + + + +/************************************************************************** + * Matrix-sparse-vector multiplication. + * + * Computes an operation of the form: + * + * y := alpha*A*x + * + * where A is a matrix and x is a sparse vector. + * + * Parameters: + * A - matrix of size n X m + * pr,ir,jc - sparse representation of the vector x, of length m + * y - output vector of length n + * alpha - real constant + * n, m - dimensions of A + * + * Note: This function re-writes the contents of y. + * + **************************************************************************/ +void mat_vec_sp(double alpha, double A[], double pr[], mwIndex ir[], mwIndex jc[], double y[], mwSize n, mwSize m); + + + +/************************************************************************** + * Matrix-transpose-sparse-vector multiplication. + * + * Computes an operation of the form: + * + * y := alpha*A'*x + * + * where A is a matrix and x is a sparse vector. + * + * Parameters: + * A - matrix of size n X m + * pr,ir,jc - sparse representation of the vector x, of length n + * y - output vector of length m + * alpha - real constant + * n, m - dimensions of A + * + * Note: This function re-writes the contents of y. + * + **************************************************************************/ +void matT_vec_sp(double alpha, double A[], double pr[], mwIndex ir[], mwIndex jc[], double y[], mwSize n, mwSize m); + + + +/************************************************************************** + * Sparse-matrix-sparse-vector multiplication. + * + * Computes an operation of the form: + * + * y := alpha*A*x + * + * where A is a sparse matrix and x is a sparse vector. + * + * Parameters: + * pr,ir,jc - sparse representation of the matrix A, of size n x m + * prx,irx,jcx - sparse representation of the vector x (of length m) + * y - output vector of length n + * alpha - real constant + * n, m - dimensions of A + * + * Note: This function re-writes the contents of y. + * + **************************************************************************/ +void mat_sp_vec_sp(double alpha, double pr[], mwIndex ir[], mwIndex jc[], double prx[], mwIndex irx[], mwIndex jcx[], double y[], mwSize n, mwSize m); + + + +/************************************************************************** + * Sparse-matrix-transpose-sparse-vector multiplication. + * + * Computes an operation of the form: + * + * y := alpha*A'*x + * + * where A is a sparse matrix and x is a sparse vector. + * + * Importnant note: this function is provided for completeness, but is NOT efficient. + * If possible, convert x to non-sparse representation and use matT_vec_sp instead. + * + * Parameters: + * pr,ir,jc - sparse representation of the matrix A, of size n x m + * prx,irx,jcx - sparse representation of the vector x (of length n) + * y - output vector of length n + * alpha - real constant + * n, m - dimensions of A + * + * Note: This function re-writes the contents of y. + * + **************************************************************************/ +void matT_sp_vec_sp(double alpha, double pr[], mwIndex ir[], mwIndex jc[], double prx[], mwIndex irx[], mwIndex jcx[], double y[], mwSize n, mwSize m); + + + +/************************************************************************** + * Matrix-matrix multiplication. + * + * Computes an operation of the form: + * + * X := alpha*A*B + * + * Parameters: + * A - matrix of size n X m + * B - matrix of size m X k + * X - output matrix of size n X k + * alpha - real constant + * n, m, k - dimensions of A, B + * + * Note: This function re-writes the contents of X. + * + **************************************************************************/ +void mat_mat(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k); + + + +/************************************************************************** + * Matrix-transpose-matrix multiplication. + * + * Computes an operation of the form: + * + * X := alpha*A*B + * + * Parameters: + * A - matrix of size n X m + * B - matrix of size m X k + * X - output matrix of size n X k + * alpha - real constant + * n, m, k - dimensions of A, B + * + * Note: This function re-writes the contents of X. + * + **************************************************************************/ +void matT_mat(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k); + + + +/************************************************************************** + * Tensor-matrix multiplication. + * + * This function accepts a 3-D tensor A of size n X m X k + * and a 2-D matrix B of size l X k. + * The function computes the 3-D tensor X of size n X m X l, where + * + * X(i,j,:) = B*A(i,j,:) + * + * for all i,j. + * + * Parameters: + * A - tensor of size n X m X k + * B - matrix of size l X k + * X - output tensor of size n X m X l + * alpha - real constant + * n, m, k, l - dimensions of A, B + * + * Note: This function re-writes the contents of X. + * + **************************************************************************/ +void tens_mat(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k, mwSize l); + + + +/************************************************************************** + * Tensor-matrix-transpose multiplication. + * + * This function accepts a 3-D tensor A of size n X m X k + * and a 2-D matrix B of size k X l. + * The function computes the 3-D tensor X of size n X m X l, where + * + * X(i,j,:) = B'*A(i,j,:) + * + * for all i,j. + * + * Parameters: + * A - tensor of size n X m X k + * B - matrix of size k X l + * X - output tensor of size n X m X l + * alpha - real constant + * n, m, k, l - dimensions of A, B + * + * Note: This function re-writes the contents of X. + * + **************************************************************************/ +void tens_matT(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k, mwSize l); + + + +/************************************************************************** + * Vector-vector sum. + * + * Computes an operation of the form: + * + * y := alpha*x + y + * + * Parameters: + * x - vector of length n + * y - output vector of length n + * alpha - real constant + * n - length of x,y + * + * Note: This function re-writes the contents of y. + * + **************************************************************************/ +void vec_sum(double alpha, double x[], double y[], mwSize n); + +/************************************************************************** + * Vector-vector scalar multiply. + * + * Computes an operation of the form: + * + * y := alpha* x.*y + * + * Parameters: + * x - vector of length n + * y - output vector of length n + * alpha - real constant + * n - length of x,y + * + * Note: This function re-writes the contents of y. + * + **************************************************************************/ + + +void vec_smult(double alpha, double x[], double y[], mwSize n); + + +/************************************************************************** + * Triangular back substitution. + * + * Solve the set of linear equations + * + * T*x = b + * + * where T is lower or upper triangular. + * + * Parameters: + * ul - 'U' for upper triangular, 'L' for lower triangular + * A - matrix of size n x m containing T + * b - vector of length k + * x - output vector of length k + * n - size of first dimension of A + * k - the size of the equation set, k<=n,m + * + * Note: + * The matrix A can be of any size n X m, as long as n,m >= k. + * Only the lower/upper triangle of the submatrix A(1:k,1:k) defines the + * matrix T (depending on the parameter ul). + * + **************************************************************************/ +void backsubst(char ul, double A[], double b[], double x[], mwSize n, mwSize k); + + + +/************************************************************************** + * Solve a set of equations using a Cholesky decomposition. + * + * Solve the set of linear equations + * + * M*x = b + * + * where M is positive definite with a known Cholesky decomposition: + * either M=L*L' (L lower triangular) or M=U'*U (U upper triangular). + * + * Parameters: + * ul - 'U' for upper triangular, 'L' for lower triangular decomposition + * A - matrix of size n x m with the Cholesky decomposition of M + * b - vector of length k + * x - output vector of length k + * n - size of first dimension of A + * k - the size of the equation set, k<=n,m + * + * Note: + * The matrix A can be of any size n X m, as long as n,m >= k. + * Only the lower/upper triangle of the submatrix A(1:k,1:k) is used as + * the Cholesky decomposition of M (depending on the parameter ul). + * + **************************************************************************/ +void cholsolve(char ul, double A[], double b[], double x[], mwSize n, mwSize k); + + + +/************************************************************************** + * Maximum absolute value. + * + * Returns the index of the coefficient with maximal absolute value in a vector. + * + * Parameters: + * x - vector of length n + * n - length of x + * + **************************************************************************/ +mwIndex maxabs(double x[], mwSize n); + + + +/************************************************************************** + * Maximum vector element. + * + * Returns the index of the maximal coefficient in a vector. + * + * Parameters: + * x - vector of length n + * n - length of x + * + **************************************************************************/ +mwIndex maxpos(double x[], mwSize n); + + + +/************************************************************************** + * Vector-vector dot product. + * + * Computes an operation of the form: + * + * c = a'*b + * + * Parameters: + * a, b - vectors of length n + * n - length of a,b + * + * Returns: The dot product c. + * + **************************************************************************/ +double dotprod(double a[], double b[], mwSize n); + + + +/************************************************************************** + * Indexed vector assignment. + * + * Perform a permutation assignment of the form + * + * y = x(ind) + * + * where ind is an array of indices to x. + * + * Parameters: + * y - output vector of length k + * x - input vector of arbitrary length + * ind - array of indices into x (indices begin at 0) + * k - length of the array ind + * + **************************************************************************/ +void vec_assign(double y[], double x[], mwIndex ind[], mwSize k); + + + +/************************************************************************** + * Matrix transpose. + * + * Computes Y := X' + * + * Parameters: + * X - input matrix of size n X m + * Y - output matrix of size m X n + * n, m - dimensions of X + * + **************************************************************************/ +void transpose(double X[], double Y[], mwSize n, mwSize m); + + + +/************************************************************************** + * Print a matrix. + * + * Parameters: + * A - matrix of size n X m + * n, m - dimensions of A + * matname - name of matrix to display + * + **************************************************************************/ +void printmat(double A[], int n, int m, char* matname); + + + +/************************************************************************** + * Print a sparse matrix. + * + * Parameters: + * A - sparse matrix of type double + * matname - name of matrix to display + * + **************************************************************************/ +void printspmat(mxArray *A, char* matname); + + +#endif +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/solvers/SMALL_ompGabor/omp2Gabor.m Mon Jul 25 17:27:05 2011 +0100 @@ -0,0 +1,223 @@ +function gamma = omp2Gabor(varargin) +%% omp2Gabor Error-constrained Orthogonal Matching Pursuit. +% +% This file is part of SMALLbox [2] and it is adaptation of Ron +% Rubenstein omp solver [1] for Gabor dictionary as defined in +% Audio Inpainting by Adler et al [3]. The dictionary is presented as +% DCT+DST and solver pics two atoms per iteration (given the frequency one +% from DCT and one from DST). For more information look in [3]. +% +% GAMMA = OMP2(D,X,G,EPSILON) solves the optimization problem +% +% min |GAMMA|_0 s.t. |X - D*GAMMA|_2 <= EPSILON +% gamma +% +% for each of the signals in X, using Batch Orthogonal Matching Pursuit. +% Here, D is a dictionary with normalized columns, X is a matrix +% containing column signals, EPSILON is the error target for each signal, +% and G is the Gramm matrix D'*D. The output GAMMA is a matrix containing +% the sparse representations as its columns. +% +% GAMMA = OMP2(D,X,[],EPSILON) performs the same operation, but without +% the matrix G, using OMP-Cholesky. This call produces the same output as +% Batch-OMP, but is significantly slower. Using this syntax is only +% recommended when available memory is too small to store G. +% +% GAMMA = OMP2(DtX,XtX,G,EPSILON) is the fastest implementation of OMP2, +% but also requires the most memory. Here, DtX stores the projections +% D'*X, and XtX is a row vector containing the squared norms of the +% signals, sum(X.*X). In this case Batch-OMP is used, but without having +% to compute D'*X and XtX in advance, which slightly improves runtime. +% Note that in general, the call +% +% GAMMA = OMP2(D'*X, sum(X.*X), G, EPSILON); +% +% will be faster than the call +% +% GAMMA = OMP2(D,X,G,EPSILON); +% +% due to optimized matrix multiplications in Matlab. However, when the +% entire matrix D'*X cannot be stored in memory, one of the other two +% versions can be used. Both compute D'*X for just one signal at a time, +% and thus require much less memory. +% +% GAMMA = OMP2(...,PARAM1,VAL1,PARAM2,VAL2,...) specifies additional +% parameters for OMP2. Available parameters are: +% +% 'gammamode' - Specifies the representation mode for GAMMA. Can be +% either 'full' or 'sparse', corresponding to a full or +% sparse matrix, respectively. By default, GAMMA is +% returned as a sparse matrix. +% 'maxatoms' - Limits the number of atoms in the representation of each +% signal. If specified, the number of atoms in each +% representation does not exceed this number, even if the +% error target is not met. Specifying maxatoms<0 implies +% no limit (default). +% 'messages' - Specifies whether progress messages should be displayed. +% When positive, this is the number of seconds between +% status prints. When negative, indicates that no messages +% should be displayed (this is the default). +% 'checkdict' - Specifies whether dictionary normalization should be +% verified. When set to 'on' (default) the dictionary +% atoms are verified to be of unit L2-norm. Setting this +% parameter to 'off' disables verification and accelerates +% function performance. Note that an unnormalized +% dictionary will produce invalid results. +% 'profile' - Can be either 'on' or 'off'. When 'on', profiling +% information is displayed at the end of the funciton +% execution. +% +% +% Summary of OMP2 versions: +% +% version | speed | memory +% ------------------------------------------------------------- +% OMP2(DtX,XtX,G,EPSILON) | very fast | very large +% OMP2(D,X,G,EPSILON) | fast | moderate +% OMP2(D,X,[],EPSILON) | very slow | small +% ------------------------------------------------------------- +% +% +% References: +% [1] M. Elad, R. Rubinstein, and M. Zibulevsky, "Efficient Implementation +% of the K-SVD Algorithm using Batch Orthogonal Matching Pursuit", +% Technical Report - CS, Technion, April 2008. +% +% See also OMP. + + +% Ron Rubinstein +% Computer Science Department +% Technion, Haifa 32000 Israel +% ronrubin@cs +% +% April 2009 + +% +% Centre for Digital Music, Queen Mary, University of London. +% This file copyright 2011 Ivan Damnjanovic. +% +% 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 2 of the +% License, or (at your option) any later version. See the file +% COPYING included with this distribution for more information. +%% + +% default options + +sparse_gamma = 1; +msgdelta = -1; +maxatoms = -1; +checkdict = 1; +profile = 0; + + +% determine number of parameters + +paramnum = 1; +while (paramnum<=nargin && ~ischar(varargin{paramnum})) + paramnum = paramnum+1; +end +paramnum = paramnum-1; + + +% parse options + +for i = paramnum+1:2:length(varargin) + paramname = varargin{i}; + paramval = varargin{i+1}; + + switch lower(paramname) + + case 'gammamode' + if (strcmpi(paramval,'sparse')) + sparse_gamma = 1; + elseif (strcmpi(paramval,'full')) + sparse_gamma = 0; + else + error('Invalid GAMMA mode'); + end + + case 'maxatoms' + maxatoms = paramval; + + case 'messages' + msgdelta = paramval; + + case 'checkdict' + if (strcmpi(paramval,'on')) + checkdict = 1; + elseif (strcmpi(paramval,'off')) + checkdict = 0; + else + error('Invalid checkdict option'); + end + + case 'profile' + if (strcmpi(paramval,'on')) + profile = 1; + elseif (strcmpi(paramval,'off')) + profile = 0; + else + error('Invalid profile mode'); + end + + otherwise + error(['Unknown option: ' paramname]); + end + +end + + +% determine call type + +if (paramnum==4) + + n1 = size(varargin{1},1); + n2 = size(varargin{2},1); + n3 = size(varargin{3},1); + + if ( (n1>1 && n2==1) || (n1==1 && n2==1 && n3==1) ) % DtX,XtX,G,EPSILON + + DtX = varargin{1}; + XtX = varargin{2}; + G = varargin{3}; + epsilon = varargin{4}; + D = []; + X = []; + + else % D,X,G,EPSILON + + D = varargin{1}; + X = varargin{2}; + G = varargin{3}; + epsilon = varargin{4}; + DtX = []; + XtX = []; + + end + +else + error('Invalid number of parameters'); +end + +G=[]; + +% verify dictionary normalization + +if (checkdict) + if (isempty(G)) + atomnorms = sum(D.*D); + else + atomnorms = diag(G); + end + if (any(abs(atomnorms-1) > 1e-2)) + error('Dictionary columns must be normalized to unit length'); + end +end + + +% omp + +gamma = omp2mexGabor(D,X,DtX,XtX,G,epsilon,sparse_gamma,msgdelta,maxatoms,profile);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/solvers/SMALL_ompGabor/omp2mex.c Mon Jul 25 17:27:05 2011 +0100 @@ -0,0 +1,156 @@ +/************************************************************************** + * + * File name: omp2mex.c + * + * Ron Rubinstein + * Computer Science Department + * Technion, Haifa 32000 Israel + * ronrubin@cs + * + * Last Updated: 18.8.2009 + * + *************************************************************************/ + +#include "ompcore.h" +#include "omputils.h" +#include "mexutils.h" + + +/* Input Arguments */ + +#define IN_D prhs[0] +#define IN_X prhs[1] +#define IN_DtX prhs[2] +#define IN_XtX prhs[3] +#define IN_G prhs[4] +#define IN_EPS prhs[5] +#define IN_SPARSE_G prhs[6] +#define IN_MSGDELTA prhs[7] +#define IN_MAXATOMS prhs[8] +#define IN_PROFILE prhs[9] + + +/* Output Arguments */ + +#define GAMMA_OUT plhs[0] + + +/***************************************************************************************/ + + +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray*prhs[]) + +{ + double *D, *x, *DtX, *XtX, *G, eps, msgdelta; + int gmode, maxatoms, profile; + mwSize m, n, L; /* D is n x m , X is n x L, DtX is m x L */ + + + /* check parameters */ + + checkmatrix(IN_D, "OMP2", "D"); + checkmatrix(IN_X, "OMP2", "X"); + checkmatrix(IN_DtX, "OMP2", "DtX"); + checkmatrix(IN_XtX, "OMP2", "XtX"); + checkmatrix(IN_G, "OMP2", "G"); + + checkscalar(IN_EPS, "OMP2", "EPSILON"); + checkscalar(IN_SPARSE_G, "OMP2", "sparse_g"); + checkscalar(IN_MSGDELTA, "OMP2", "msgdelta"); + checkscalar(IN_MAXATOMS, "OMP2", "maxatoms"); + checkscalar(IN_PROFILE, "OMP2", "profile"); + + + /* get parameters */ + + x = D = DtX = XtX = G = 0; + + if (!mxIsEmpty(IN_D)) + D = mxGetPr(IN_D); + + if (!mxIsEmpty(IN_X)) + x = mxGetPr(IN_X); + + if (!mxIsEmpty(IN_DtX)) + DtX = mxGetPr(IN_DtX); + + if (!mxIsEmpty(IN_XtX)) + XtX = mxGetPr(IN_XtX); + + if (!mxIsEmpty(IN_G)) + G = mxGetPr(IN_G); + + eps = mxGetScalar(IN_EPS); + if ((int)(mxGetScalar(IN_SPARSE_G)+1e-2)) { + gmode = SPARSE_GAMMA; + } + else { + gmode = FULL_GAMMA; + } + msgdelta = mxGetScalar(IN_MSGDELTA); + if (mxGetScalar(IN_MAXATOMS) < -1e-5) { + maxatoms = -1; + } + else { + maxatoms = (int)(mxGetScalar(IN_MAXATOMS)+1e-2); + } + profile = (int)(mxGetScalar(IN_PROFILE)+1e-2); + + + /* check sizes */ + + if (D && x) { + n = mxGetM(IN_D); + m = mxGetN(IN_D); + L = mxGetN(IN_X); + + if (mxGetM(IN_X) != n) { + mexErrMsgTxt("D and X have incompatible sizes."); + } + + if (G) { + if (mxGetN(IN_G)!=mxGetM(IN_G)) { + mexErrMsgTxt("G must be a square matrix."); + } + if (mxGetN(IN_G) != m) { + mexErrMsgTxt("D and G have incompatible sizes."); + } + } + } + + else if (DtX && XtX) { + m = mxGetM(IN_DtX); + L = mxGetN(IN_DtX); + + /* set n to an arbitrary value that is at least the max possible number of selected atoms */ + + if (maxatoms>0) { + n = maxatoms; + } + else { + n = m; + } + + if ( !(mxGetM(IN_XtX)==L && mxGetN(IN_XtX)==1) && !(mxGetM(IN_XtX)==1 && mxGetN(IN_XtX)==L) ) { + mexErrMsgTxt("DtX and XtX have incompatible sizes."); + } + + if (mxGetN(IN_G)!=mxGetM(IN_G)) { + mexErrMsgTxt("G must be a square matrix."); + } + if (mxGetN(IN_G) != m) { + mexErrMsgTxt("DtX and G have incompatible sizes."); + } + } + + else { + mexErrMsgTxt("Either D and X, or DtX and XtX, must be specified."); + } + + + /* Do OMP! */ + + GAMMA_OUT = ompcore(D, x, DtX, XtX, G, n, m, L, maxatoms, eps, gmode, profile, msgdelta, 1); + + return; +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/solvers/SMALL_ompGabor/omp2mex.m Mon Jul 25 17:27:05 2011 +0100 @@ -0,0 +1,23 @@ +%This is the Matlab interface to the OMP2 MEX implementation. +%The function is not for independent use, only through omp2.m. + + +%OMP2MEX Matlab interface to the OMP2 MEX implementation. +% GAMMA = OMP2MEX(D,X,DtX,XtX,G,EPSILON,SPARSE_G,MSGDELTA,MAXATOMS,PROFILE) +% invokes the OMP2 MEX function according to the specified parameters. Not +% all the parameters are required. Those among D, X, DtX, XtX and G which +% are not specified should be passed as []. +% +% EPSILON - the target error. +% SPARSE_G - returns a sparse GAMMA when nonzero, full GAMMA when zero. +% MSGDELTA - the delay in secs between messages. Zero means no messages. +% MAXATOMS - the max number of atoms per signal, negative for no max. +% PROFILE - nonzero means that profiling information should be printed. + + +% Ron Rubinstein +% Computer Science Department +% Technion, Haifa 32000 Israel +% ronrubin@cs +% +% April 2009
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/solvers/SMALL_ompGabor/omp2mexGabor.c Mon Jul 25 17:27:05 2011 +0100 @@ -0,0 +1,156 @@ +/************************************************************************** + * + * File name: omp2mex.c + * + * Ron Rubinstein + * Computer Science Department + * Technion, Haifa 32000 Israel + * ronrubin@cs + * + * Last Updated: 18.8.2009 + * + *************************************************************************/ + +#include "ompcoreGabor.h" +#include "omputils.h" +#include "mexutils.h" + + +/* Input Arguments */ + +#define IN_D prhs[0] +#define IN_X prhs[1] +#define IN_DtX prhs[2] +#define IN_XtX prhs[3] +#define IN_G prhs[4] +#define IN_EPS prhs[5] +#define IN_SPARSE_G prhs[6] +#define IN_MSGDELTA prhs[7] +#define IN_MAXATOMS prhs[8] +#define IN_PROFILE prhs[9] + + +/* Output Arguments */ + +#define GAMMA_OUT plhs[0] + + +/***************************************************************************************/ + + +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray*prhs[]) + +{ + double *D, *x, *DtX, *XtX, *G, eps, msgdelta; + int gmode, maxatoms, profile; + mwSize m, n, L; /* D is n x m , X is n x L, DtX is m x L */ + + + /* check parameters */ + + checkmatrix(IN_D, "OMP2", "D"); + checkmatrix(IN_X, "OMP2", "X"); + checkmatrix(IN_DtX, "OMP2", "DtX"); + checkmatrix(IN_XtX, "OMP2", "XtX"); + checkmatrix(IN_G, "OMP2", "G"); + + checkscalar(IN_EPS, "OMP2", "EPSILON"); + checkscalar(IN_SPARSE_G, "OMP2", "sparse_g"); + checkscalar(IN_MSGDELTA, "OMP2", "msgdelta"); + checkscalar(IN_MAXATOMS, "OMP2", "maxatoms"); + checkscalar(IN_PROFILE, "OMP2", "profile"); + + + /* get parameters */ + + x = D = DtX = XtX = G = 0; + + if (!mxIsEmpty(IN_D)) + D = mxGetPr(IN_D); + + if (!mxIsEmpty(IN_X)) + x = mxGetPr(IN_X); + + if (!mxIsEmpty(IN_DtX)) + DtX = mxGetPr(IN_DtX); + + if (!mxIsEmpty(IN_XtX)) + XtX = mxGetPr(IN_XtX); + + if (!mxIsEmpty(IN_G)) + G = mxGetPr(IN_G); + + eps = mxGetScalar(IN_EPS); + if ((int)(mxGetScalar(IN_SPARSE_G)+1e-2)) { + gmode = SPARSE_GAMMA; + } + else { + gmode = FULL_GAMMA; + } + msgdelta = mxGetScalar(IN_MSGDELTA); + if (mxGetScalar(IN_MAXATOMS) < -1e-5) { + maxatoms = -1; + } + else { + maxatoms = (int)(mxGetScalar(IN_MAXATOMS)+1e-2); + } + profile = (int)(mxGetScalar(IN_PROFILE)+1e-2); + + + /* check sizes */ + + if (D && x) { + n = mxGetM(IN_D); + m = mxGetN(IN_D); + L = mxGetN(IN_X); + + if (mxGetM(IN_X) != n) { + mexErrMsgTxt("D and X have incompatible sizes."); + } + + if (G) { + if (mxGetN(IN_G)!=mxGetM(IN_G)) { + mexErrMsgTxt("G must be a square matrix."); + } + if (mxGetN(IN_G) != m) { + mexErrMsgTxt("D and G have incompatible sizes."); + } + } + } + + else if (DtX && XtX) { + m = mxGetM(IN_DtX); + L = mxGetN(IN_DtX); + + /* set n to an arbitrary value that is at least the max possible number of selected atoms */ + + if (maxatoms>0) { + n = maxatoms; + } + else { + n = m; + } + + if ( !(mxGetM(IN_XtX)==L && mxGetN(IN_XtX)==1) && !(mxGetM(IN_XtX)==1 && mxGetN(IN_XtX)==L) ) { + mexErrMsgTxt("DtX and XtX have incompatible sizes."); + } + + if (mxGetN(IN_G)!=mxGetM(IN_G)) { + mexErrMsgTxt("G must be a square matrix."); + } + if (mxGetN(IN_G) != m) { + mexErrMsgTxt("DtX and G have incompatible sizes."); + } + } + + else { + mexErrMsgTxt("Either D and X, or DtX and XtX, must be specified."); + } + + + /* Do OMP! */ + + GAMMA_OUT = ompcoreGabor(D, x, DtX, XtX, G, n, m, L, maxatoms, eps, gmode, profile, msgdelta, 1); + + return; +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/solvers/SMALL_ompGabor/omp2mexGabor.m Mon Jul 25 17:27:05 2011 +0100 @@ -0,0 +1,23 @@ +%This is the Matlab interface to the OMP2 MEX implementation. +%The function is not for independent use, only through omp2.m. + + +%OMP2MEX Matlab interface to the OMP2 MEX implementation. +% GAMMA = OMP2MEXGABOR(D,X,DtX,XtX,G,EPSILON,SPARSE_G,MSGDELTA,MAXATOMS,PROFILE) +% invokes the OMP2 MEX function according to the specified parameters. Not +% all the parameters are required. Those among D, X, DtX, XtX and G which +% are not specified should be passed as []. +% +% EPSILON - the target error. +% SPARSE_G - returns a sparse GAMMA when nonzero, full GAMMA when zero. +% MSGDELTA - the delay in secs between messages. Zero means no messages. +% MAXATOMS - the max number of atoms per signal, negative for no max. +% PROFILE - nonzero means that profiling information should be printed. + + +% Ron Rubinstein +% Computer Science Department +% Technion, Haifa 32000 Israel +% ronrubin@cs +% +% April 2009
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/solvers/SMALL_ompGabor/ompGabor.m Mon Jul 25 17:27:05 2011 +0100 @@ -0,0 +1,180 @@ +function gamma = omp(varargin) +%OMP Sparsity-constrained Orthogonal Matching Pursuit. +% GAMMA = OMP(D,X,G,T) solves the optimization problem +% +% min |X - D*GAMMA|_2 s.t. |GAMMA|_0 <= T +% gamma +% +% for each of the signals in X, using Batch Orthogonal Matching Pursuit. +% Here, D is a dictionary with normalized columns, X is a matrix +% containing column signals, T is the # of non-zeros in each signal +% representation, and G is the Gramm matrix D'*D. The output GAMMA is a +% matrix containing the sparse representations as its columns. +% +% GAMMA = OMP(D,X,[],T) performs the same operation, but without the +% matrix G, using OMP-Cholesky. This call produces the same output as +% Batch-OMP, but is significantly slower. Using this syntax is only +% recommended when available memory is too small to store G. +% +% GAMMA = OMP(DtX,G,T) is the fastest implementation of OMP, but also +% requires the most memory. Here, DtX stores the projections D'*X. In this +% case Batch-OMP is used, but without having to compute D'*X in advance, +% which slightly improves runtime. Note that in general, the call +% +% GAMMA = OMP(D'*X,G,T); +% +% will be faster than the call +% +% GAMMA = OMP(D,X,G,T); +% +% due to optimized matrix multiplications in Matlab. However, when the +% entire matrix D'*X cannot be stored in memory, one of the other two +% versions can be used. Both compute D'*X for just one signal at a time, +% and thus require much less memory. +% +% GAMMA = OMP(...,PARAM1,VAL1,PARAM2,VAL2,...) specifies additional +% parameters for OMP. Available parameters are: +% +% 'gammamode' - Specifies the representation mode for GAMMA. Can be +% either 'full' or 'sparse', corresponding to a full or +% sparse matrix, respectively. By default, GAMMA is +% returned as a sparse matrix. +% 'messages' - Specifies whether progress messages should be displayed. +% When positive, this is the number of seconds between +% status prints. When negative, indicates that no messages +% should be displayed (this is the default). +% 'checkdict' - Specifies whether dictionary normalization should be +% verified. When set to 'on' (default) the dictionary +% atoms are verified to be of unit L2-norm. Setting this +% parameter to 'off' disables verification and accelerates +% function performance. Note that an unnormalized +% dictionary will produce invalid results. +% 'profile' - Can be either 'on' or 'off'. When 'on', profiling +% information is displayed at the end of the funciton +% execution. +% +% +% Summary of OMP versions: +% +% version | speed | memory +% -------------------------------------------------- +% OMP(DtX,G,T) | very fast | very large +% OMP(D,X,G,T) | fast | moderate +% OMP(D,X,[],T) | very slow | small +% -------------------------------------------------- +% +% +% References: +% [1] M. Elad, R. Rubinstein, and M. Zibulevsky, "Efficient Implementation +% of the K-SVD Algorithm using Batch Orthogonal Matching Pursuit", +% Technical Report - CS, Technion, April 2008. +% +% See also OMP2. + + +% Ron Rubinstein +% Computer Science Department +% Technion, Haifa 32000 Israel +% ronrubin@cs +% +% April 2009 + + +% default options + +sparse_gamma = 1; +msgdelta = -1; +checkdict = 1; +profile = 0; + + +% determine number of parameters + +paramnum = 1; +while (paramnum<=nargin && ~ischar(varargin{paramnum})) + paramnum = paramnum+1; +end +paramnum = paramnum-1; + + +% parse options + +for i = paramnum+1:2:length(varargin) + paramname = varargin{i}; + paramval = varargin{i+1}; + + switch lower(paramname) + + case 'gammamode' + if (strcmpi(paramval,'sparse')) + sparse_gamma = 1; + elseif (strcmpi(paramval,'full')) + sparse_gamma = 0; + else + error('Invalid GAMMA mode'); + end + + case 'messages' + msgdelta = paramval; + + case 'checkdict' + if (strcmpi(paramval,'on')) + checkdict = 1; + elseif (strcmpi(paramval,'off')) + checkdict = 0; + else + error('Invalid checkdict option'); + end + + case 'profile' + if (strcmpi(paramval,'on')) + profile = 1; + elseif (strcmpi(paramval,'off')) + profile = 0; + else + error('Invalid profile mode'); + end + + otherwise + error(['Unknown option: ' paramname]); + end + +end + + +% determine call type + +if (paramnum==3) + DtX = varargin{1}; + G = varargin{2}; + T = varargin{3}; + D = []; + X = []; +elseif (paramnum==4) + D = varargin{1}; + X = varargin{2}; + G = varargin{3}; + T = varargin{4}; + DtX = []; +else + error('Invalid number of parameters'); +end + + +% verify dictionary normalization + +if (checkdict) + if (isempty(G)) + atomnorms = sum(D.*D); + else + atomnorms = diag(G); + end + if (any(abs(atomnorms-1) > 1e-2)) + error('Dictionary columns must be normalized to unit length'); + end +end + + +% omp + +gamma = ompmexGabor(D,X,DtX,G,T,sparse_gamma,msgdelta,profile);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/solvers/SMALL_ompGabor/ompcore.c Mon Jul 25 17:27:05 2011 +0100 @@ -0,0 +1,409 @@ +/************************************************************************** + * + * File name: ompcore.c + * + * Ron Rubinstein + * Computer Science Department + * Technion, Haifa 32000 Israel + * ronrubin@cs + * + * Last Updated: 25.8.2009 + * + *************************************************************************/ + + +#include "ompcore.h" +#include "omputils.h" +#include "ompprof.h" +#include "myblas.h" +#include <math.h> +#include <string.h> + + + +/****************************************************************************** + * * + * Batch-OMP Implementation * + * * + ******************************************************************************/ + +mxArray* ompcore(double D[], double x[], double DtX[], double XtX[], double G[], mwSize n, mwSize m, mwSize L, + int T, double eps, int gamma_mode, int profile, double msg_delta, int erroromp) +{ + + profdata pd; + mxArray *Gamma; + mwIndex i, j, signum, pos, *ind, *gammaIr, *gammaJc, gamma_count; + mwSize allocated_coefs, allocated_cols; + int DtX_specified, XtX_specified, batchomp, standardomp, *selected_atoms; + double *alpha, *r, *Lchol, *c, *Gsub, *Dsub, sum, *gammaPr, *tempvec1, *tempvec2; + double eps2, resnorm, delta, deltaprev, secs_remain; + int mins_remain, hrs_remain; + clock_t lastprint_time, starttime; + + + + /*** status flags ***/ + + DtX_specified = (DtX!=0); /* indicates whether D'*x was provided */ + XtX_specified = (XtX!=0); /* indicates whether sum(x.*x) was provided */ + + standardomp = (G==0); /* batch-omp or standard omp are selected depending on availability of G */ + batchomp = !standardomp; + + + + /*** allocate output matrix ***/ + + + if (gamma_mode == FULL_GAMMA) { + + /* allocate full matrix of size m X L */ + + Gamma = mxCreateDoubleMatrix(m, L, mxREAL); + gammaPr = mxGetPr(Gamma); + gammaIr = 0; + gammaJc = 0; + } + else { + + /* allocate sparse matrix with room for allocated_coefs nonzeros */ + + /* for error-omp, begin with L*sqrt(n)/2 allocated nonzeros, otherwise allocate L*T nonzeros */ + allocated_coefs = erroromp ? (mwSize)(ceil(L*sqrt((double)n)/2.0) + 1.01) : L*T; + Gamma = mxCreateSparse(m, L, allocated_coefs, mxREAL); + gammaPr = mxGetPr(Gamma); + gammaIr = mxGetIr(Gamma); + gammaJc = mxGetJc(Gamma); + gamma_count = 0; + gammaJc[0] = 0; + } + + + /*** helper arrays ***/ + + alpha = (double*)mxMalloc(m*sizeof(double)); /* contains D'*residual */ + ind = (mwIndex*)mxMalloc(n*sizeof(mwIndex)); /* indices of selected atoms */ + selected_atoms = (int*)mxMalloc(m*sizeof(int)); /* binary array with 1's for selected atoms */ + c = (double*)mxMalloc(n*sizeof(double)); /* orthogonal projection result */ + + /* current number of columns in Dsub / Gsub / Lchol */ + allocated_cols = erroromp ? (mwSize)(ceil(sqrt((double)n)/2.0) + 1.01) : T; + + /* Cholesky decomposition of D_I'*D_I */ + Lchol = (double*)mxMalloc(n*allocated_cols*sizeof(double)); + + /* temporary vectors for various computations */ + tempvec1 = (double*)mxMalloc(m*sizeof(double)); + tempvec2 = (double*)mxMalloc(m*sizeof(double)); + + if (batchomp) { + /* matrix containing G(:,ind) - the columns of G corresponding to the selected atoms, in order of selection */ + Gsub = (double*)mxMalloc(m*allocated_cols*sizeof(double)); + } + else { + /* matrix containing D(:,ind) - the selected atoms from D, in order of selection */ + Dsub = (double*)mxMalloc(n*allocated_cols*sizeof(double)); + + /* stores the residual */ + r = (double*)mxMalloc(n*sizeof(double)); + } + + if (!DtX_specified) { + /* contains D'*x for the current signal */ + DtX = (double*)mxMalloc(m*sizeof(double)); + } + + + + /*** initializations for error omp ***/ + + if (erroromp) { + eps2 = eps*eps; /* compute eps^2 */ + if (T<0 || T>n) { /* unspecified max atom num - set max atoms to n */ + T = n; + } + } + + + + /*** initialize timers ***/ + + initprofdata(&pd); /* initialize profiling counters */ + starttime = clock(); /* record starting time for eta computations */ + lastprint_time = starttime; /* time of last status display */ + + + + /********************** perform omp for each signal **********************/ + + + + for (signum=0; signum<L; ++signum) { + + + /* initialize residual norm and deltaprev for error-omp */ + + if (erroromp) { + if (XtX_specified) { + resnorm = XtX[signum]; + } + else { + resnorm = dotprod(x+n*signum, x+n*signum, n); + addproftime(&pd, XtX_TIME); + } + deltaprev = 0; /* delta tracks the value of gamma'*G*gamma */ + } + else { + /* ignore residual norm stopping criterion */ + eps2 = 0; + resnorm = 1; + } + + + if (resnorm>eps2 && T>0) { + + /* compute DtX */ + + if (!DtX_specified) { + matT_vec(1, D, x+n*signum, DtX, n, m); + addproftime(&pd, DtX_TIME); + } + + + /* initialize alpha := DtX */ + + memcpy(alpha, DtX + m*signum*DtX_specified, m*sizeof(double)); + + + /* mark all atoms as unselected */ + + for (i=0; i<m; ++i) { + selected_atoms[i] = 0; + } + + } + + + /* main loop */ + + i=0; + while (resnorm>eps2 && i<T) { + + /* index of next atom */ + + pos = maxabs(alpha, m); + addproftime(&pd, MAXABS_TIME); + + + /* stop criterion: selected same atom twice, or inner product too small */ + + if (selected_atoms[pos] || alpha[pos]*alpha[pos]<1e-14) { + break; + } + + + /* mark selected atom */ + + ind[i] = pos; + selected_atoms[pos] = 1; + + + /* matrix reallocation */ + + if (erroromp && i>=allocated_cols) { + + allocated_cols = (mwSize)(ceil(allocated_cols*MAT_INC_FACTOR) + 1.01); + + Lchol = (double*)mxRealloc(Lchol,n*allocated_cols*sizeof(double)); + + batchomp ? (Gsub = (double*)mxRealloc(Gsub,m*allocated_cols*sizeof(double))) : + (Dsub = (double*)mxRealloc(Dsub,n*allocated_cols*sizeof(double))) ; + } + + + /* append column to Gsub or Dsub */ + + if (batchomp) { + memcpy(Gsub+i*m, G+pos*m, m*sizeof(double)); + } + else { + memcpy(Dsub+i*n, D+pos*n, n*sizeof(double)); + } + + + /*** Cholesky update ***/ + + if (i==0) { + *Lchol = 1; + } + else { + + /* incremental Cholesky decomposition: compute next row of Lchol */ + + if (standardomp) { + matT_vec(1, Dsub, D+n*pos, tempvec1, n, i); /* compute tempvec1 := Dsub'*d where d is new atom */ + addproftime(&pd, DtD_TIME); + } + else { + vec_assign(tempvec1, Gsub+i*m, ind, i); /* extract tempvec1 := Gsub(ind,i) */ + } + backsubst('L', Lchol, tempvec1, tempvec2, n, i); /* compute tempvec2 = Lchol \ tempvec1 */ + for (j=0; j<i; ++j) { /* write tempvec2 to end of Lchol */ + Lchol[j*n+i] = tempvec2[j]; + } + + /* compute Lchol(i,i) */ + sum = 0; + for (j=0; j<i; ++j) { /* compute sum of squares of last row without Lchol(i,i) */ + sum += SQR(Lchol[j*n+i]); + } + if ( (1-sum) <= 1e-14 ) { /* Lchol(i,i) is zero => selected atoms are dependent */ + break; + } + Lchol[i*n+i] = sqrt(1-sum); + } + + addproftime(&pd, LCHOL_TIME); + + i++; + + + /* perform orthogonal projection and compute sparse coefficients */ + + vec_assign(tempvec1, DtX + m*signum*DtX_specified, ind, i); /* extract tempvec1 = DtX(ind) */ + cholsolve('L', Lchol, tempvec1, c, n, i); /* solve LL'c = tempvec1 for c */ + addproftime(&pd, COMPCOEF_TIME); + + + /* update alpha = D'*residual */ + + if (standardomp) { + mat_vec(-1, Dsub, c, r, n, i); /* compute r := -Dsub*c */ + vec_sum(1, x+n*signum, r, n); /* compute r := x+r */ + + + /*memcpy(r, x+n*signum, n*sizeof(double)); /* assign r := x */ + /*mat_vec1(-1, Dsub, c, 1, r, n, i); /* compute r := r-Dsub*c */ + + addproftime(&pd, COMPRES_TIME); + matT_vec(1, D, r, alpha, n, m); /* compute alpha := D'*r */ + addproftime(&pd, DtR_TIME); + + /* update residual norm */ + if (erroromp) { + resnorm = dotprod(r, r, n); + addproftime(&pd, UPDATE_RESNORM_TIME); + } + } + else { + mat_vec(1, Gsub, c, tempvec1, m, i); /* compute tempvec1 := Gsub*c */ + memcpy(alpha, DtX + m*signum*DtX_specified, m*sizeof(double)); /* set alpha = D'*x */ + vec_sum(-1, tempvec1, alpha, m); /* compute alpha := alpha - tempvec1 */ + addproftime(&pd, UPDATE_DtR_TIME); + + /* update residual norm */ + if (erroromp) { + vec_assign(tempvec2, tempvec1, ind, i); /* assign tempvec2 := tempvec1(ind) */ + delta = dotprod(c,tempvec2,i); /* compute c'*tempvec2 */ + resnorm = resnorm - delta + deltaprev; /* residual norm update */ + deltaprev = delta; + addproftime(&pd, UPDATE_RESNORM_TIME); + } + } + } + + + /*** generate output vector gamma ***/ + + if (gamma_mode == FULL_GAMMA) { /* write the coefs in c to their correct positions in gamma */ + for (j=0; j<i; ++j) { + gammaPr[m*signum + ind[j]] = c[j]; + } + } + else { + /* sort the coefs by index before writing them to gamma */ + quicksort(ind,c,i); + addproftime(&pd, INDEXSORT_TIME); + + /* gamma is full - reallocate */ + if (gamma_count+i >= allocated_coefs) { + + while(gamma_count+i >= allocated_coefs) { + allocated_coefs = (mwSize)(ceil(GAMMA_INC_FACTOR*allocated_coefs) + 1.01); + } + + mxSetNzmax(Gamma, allocated_coefs); + mxSetPr(Gamma, mxRealloc(gammaPr, allocated_coefs*sizeof(double))); + mxSetIr(Gamma, mxRealloc(gammaIr, allocated_coefs*sizeof(mwIndex))); + + gammaPr = mxGetPr(Gamma); + gammaIr = mxGetIr(Gamma); + } + + /* append coefs to gamma and update the indices */ + for (j=0; j<i; ++j) { + gammaPr[gamma_count] = c[j]; + gammaIr[gamma_count] = ind[j]; + gamma_count++; + } + gammaJc[signum+1] = gammaJc[signum] + i; + } + + + + /*** display status messages ***/ + + if (msg_delta>0 && (clock()-lastprint_time)/(double)CLOCKS_PER_SEC >= msg_delta) + { + lastprint_time = clock(); + + /* estimated remainig time */ + secs2hms( ((L-signum-1)/(double)(signum+1)) * ((lastprint_time-starttime)/(double)CLOCKS_PER_SEC) , + &hrs_remain, &mins_remain, &secs_remain); + + mexPrintf("omp: signal %d / %d, estimated remaining time: %02d:%02d:%05.2f\n", + signum+1, L, hrs_remain, mins_remain, secs_remain); + mexEvalString("drawnow;"); + } + + } + + /* end omp */ + + + + /*** print final messages ***/ + + if (msg_delta>0) { + mexPrintf("omp: signal %d / %d\n", signum, L); + } + + if (profile) { + printprofinfo(&pd, erroromp, batchomp, L); + } + + + + /* free memory */ + + if (!DtX_specified) { + mxFree(DtX); + } + if (standardomp) { + mxFree(r); + mxFree(Dsub); + } + else { + mxFree(Gsub); + } + mxFree(tempvec2); + mxFree(tempvec1); + mxFree(Lchol); + mxFree(c); + mxFree(selected_atoms); + mxFree(ind); + mxFree(alpha); + + return Gamma; +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/solvers/SMALL_ompGabor/ompcore.h Mon Jul 25 17:27:05 2011 +0100 @@ -0,0 +1,80 @@ +/************************************************************************** + * + * File name: ompcore.h + * + * Ron Rubinstein + * Computer Science Department + * Technion, Haifa 32000 Israel + * ronrubin@cs + * + * Last Updated: 18.8.2009 + * + * Contains the core implementation of Batch-OMP / OMP-Cholesky. + * + *************************************************************************/ + + +#ifndef __OMP_CORE_H__ +#define __OMP_CORE_H__ + + +#include "mex.h" + + + +/************************************************************************** + * Perform Batch-OMP or OMP-Cholesky on a specified set of signals, using + * either a fixed number of atoms or an error bound. + * + * Parameters (not all required): + * + * D - the dictionary, of size n X m + * x - the signals, of size n X L + * DtX - D'*x, of size m X L + * XtX - squared norms of the signals in x, sum(x.*x), of length L + * G - D'*D, of size m X m + * T - target sparsity, or maximal number of atoms for error-based OMP + * eps - target residual norm for error-based OMP + * gamma_mode - one of the constants FULL_GAMMA or SPARSE_GAMMA + * profile - if non-zero, profiling info is printed + * msg_delta - positive: the # of seconds between status prints, otherwise: nothing is printed + * erroromp - if nonzero indicates error-based OMP, otherwise fixed sparsity OMP + * + * Usage: + * + * The function can be called using different parameters, and will have + * different complexity depending on the parameters specified. Arrays which + * are not specified should be passed as null (0). When G is specified, + * Batch-OMP is performed. Otherwise, OMP-Cholesky is performed. + * + * Fixed-sparsity usage: + * --------------------- + * Either DtX, or D and x, must be specified. Specifying DtX is more efficient. + * XtX does not need to be specified. + * When D and x are specified, G is not required. However, not providing G + * will significantly degrade efficiency. + * The number of atoms must be specified in T. The value of eps is ignored. + * Finally, set erroromp to 0. + * + * Error-OMP usage: + * ---------------- + * Either DtX and Xtx, or D and x, must be specified. Specifying DtX and XtX + * is more efficient. + * When D and x are specified, G is not required. However, not providing G + * will significantly degrade efficiency. + * The target error must be specified in eps. A hard limit on the number + * of atoms can also be specified via the parameter T. Otherwise, T should + * be negative. Finally, set erroromp to nonzero. + * + * + * Returns: + * An mxArray containing the sparse representations of the signals in x + * (allocated using the appropriate mxCreateXXX() function). + * The array is either full or sparse, depending on gamma_mode. + * + **************************************************************************/ +mxArray* ompcore(double D[], double x[], double DtX[], double XtX[], double G[], mwSize n, mwSize m, mwSize L, + int T, double eps, int gamma_mode, int profile, double msg_delta, int erroromp); + + +#endif
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/solvers/SMALL_ompGabor/ompcoreGabor.c Mon Jul 25 17:27:05 2011 +0100 @@ -0,0 +1,465 @@ +/************************************************************************** + * + * File name: ompcoreGabor.c + * + * Ron Rubinstein + * Computer Science Department + * Technion, Haifa 32000 Israel + * ronrubin@cs + * + * Last Updated: 25.8.2009 + * + * Modified by Ivan damnjanovic July 2011 + * Takes to atoms per iteration. It should be used for Gabor dictionaries + * as specified in + * "Audio Inpainting" Amir Adler, Valentin Emiya, Maria G. Jafari, + * Michael Elad, Remi Gribonval and Mark D. Plumbley + * Draft version: March 6, 2011 + * + *************************************************************************/ + + +#include "ompcoreGabor.h" +#include "omputils.h" +#include "ompprof.h" +#include "myblas.h" +#include <math.h> +#include <string.h> + + + +/****************************************************************************** + * * + * Batch-OMP Implementation * + * * + ******************************************************************************/ + +mxArray* ompcoreGabor(double D[], double x[], double DtX[], double XtX[], double G[], mwSize n, mwSize m, mwSize L, + int T, double eps, int gamma_mode, int profile, double msg_delta, int erroromp) +{ + + profdata pd; + mxArray *Gamma; + mwIndex i, j, k, signum, pos, *ind, *gammaIr, *gammaJc, gamma_count; + mwSize allocated_coefs, allocated_cols; + int DtX_specified, XtX_specified, batchomp, standardomp, *selected_atoms; + double *proj, *proj1, *proj2, *D1, *D2, *D1D2, *n12, *alpha, *beta, *error; + double *r, *Lchol, *c, *Gsub, *Dsub, sum, *gammaPr, *tempvec1, *tempvec2; + double eps2, resnorm, delta, deltaprev, secs_remain; + int mins_remain, hrs_remain; + clock_t lastprint_time, starttime; + + + + /*** status flags ***/ + + DtX_specified = (DtX!=0); /* indicates whether D'*x was provided */ + XtX_specified = (XtX!=0); /* indicates whether sum(x.*x) was provided */ + + standardomp = (G==0); /* batch-omp or standard omp are selected depending on availability of G */ + batchomp = !standardomp; + + + + /*** allocate output matrix ***/ + + + if (gamma_mode == FULL_GAMMA) { + + /* allocate full matrix of size m X L */ + + Gamma = mxCreateDoubleMatrix(m, L, mxREAL); + gammaPr = mxGetPr(Gamma); + gammaIr = 0; + gammaJc = 0; + } + else { + + /* allocate sparse matrix with room for allocated_coefs nonzeros */ + + /* for error-omp, begin with L*sqrt(n)/2 allocated nonzeros, otherwise allocate L*T nonzeros */ + allocated_coefs = erroromp ? (mwSize)(ceil(L*sqrt((double)n)/2.0) + 1.01) : L*T; + Gamma = mxCreateSparse(m, L, allocated_coefs, mxREAL); + gammaPr = mxGetPr(Gamma); + gammaIr = mxGetIr(Gamma); + gammaJc = mxGetJc(Gamma); + gamma_count = 0; + gammaJc[0] = 0; + } + + + /*** helper arrays ***/ + /* Ivan Damnjanovic July 2011*/ + proj = (double*)mxMalloc(m*sizeof(double)); + proj1 = (double*)mxMalloc(m/2*sizeof(double)); + proj2 = (double*)mxMalloc(m/2*sizeof(double)); + D1 = (double*)mxMalloc(n*m/2*sizeof(double)); + D2 = (double*)mxMalloc(n*m/2*sizeof(double)); + memcpy(D1, D , n*m/2*sizeof(double)); + memcpy(D2, D+n*m/2, n*m/2*sizeof(double)); + + D1D2 = (double*)mxMalloc(m/2*sizeof(double)); + n12 = (double*)mxMalloc(m/2*sizeof(double)); + + vec_smult(1,D2, D1, n*m/2); + + for (i=0; i<m/2; i++) { + D1D2[i]=0; + n12[i]=0; + for (j=0; j<n; j++) { + D1D2[i] += D1[i*n+j]; + } + n12[i]=1/(1-D1D2[i]*D1D2[i]); + } + + memcpy(D1, D , n*m/2*sizeof(double)); + + alpha = (double*)mxMalloc(m/2*sizeof(double)); /* contains D'*residual */ + beta = (double*)mxMalloc(m/2*sizeof(double)); + error = (double*)mxMalloc(m/2*sizeof(double)); + + ind = (mwIndex*)mxMalloc(m*sizeof(mwIndex)); /* indices of selected atoms */ + selected_atoms = (int*)mxMalloc(m*sizeof(int)); /* binary array with 1's for selected atoms */ + c = (double*)mxMalloc(n*sizeof(double)); /* orthogonal projection result */ + + /* current number of columns in Dsub / Gsub / Lchol */ + allocated_cols = erroromp ? (mwSize)(ceil(sqrt((double)n)/2.0) + 1.01) : T; + + /* Cholesky decomposition of D_I'*D_I */ + Lchol = (double*)mxMalloc(n*allocated_cols*sizeof(double)); + + /* temporary vectors for various computations */ + tempvec1 = (double*)mxMalloc(m*sizeof(double)); + tempvec2 = (double*)mxMalloc(m*sizeof(double)); + + if (batchomp) { + /* matrix containing G(:,ind) - the columns of G corresponding to the selected atoms, in order of selection */ + Gsub = (double*)mxMalloc(m*allocated_cols*sizeof(double)); + } + else { + /* matrix containing D(:,ind) - the selected atoms from D, in order of selection */ + Dsub = (double*)mxMalloc(n*allocated_cols*sizeof(double)); + + /* stores the residual */ + r = (double*)mxMalloc(n*sizeof(double)); + } + + if (!DtX_specified) { + /* contains D'*x for the current signal */ + DtX = (double*)mxMalloc(m*sizeof(double)); + } + + + + /*** initializations for error omp ***/ + + if (erroromp) { + eps2 = eps*eps; /* compute eps^2 */ + if (T<0 || T>n) { /* unspecified max atom num - set max atoms to n */ + T = n; + } + } + + + + /*** initialize timers ***/ + + initprofdata(&pd); /* initialize profiling counters */ + starttime = clock(); /* record starting time for eta computations */ + lastprint_time = starttime; /* time of last status display */ + + + + /********************** perform omp for each signal **********************/ + + + + for (signum=0; signum<L; ++signum) { + + + /* initialize residual norm and deltaprev for error-omp */ + + if (erroromp) { + if (XtX_specified) { + resnorm = XtX[signum]; + } + else { + resnorm = dotprod(x+n*signum, x+n*signum, n); + addproftime(&pd, XtX_TIME); + } + deltaprev = 0; /* delta tracks the value of gamma'*G*gamma */ + } + else { + /* ignore residual norm stopping criterion */ + eps2 = 0; + resnorm = 1; + } + + + if (resnorm>eps2 && T>0) { + + /* compute DtX */ + + if (!DtX_specified) { + matT_vec(1, D, x+n*signum, DtX, n, m); + addproftime(&pd, DtX_TIME); + memcpy(r , x+n*signum, n*sizeof(double)); + } + + + /* initialize projections to D1 and D2 := DtX */ + + memcpy(proj, DtX + m*signum*DtX_specified, m*sizeof(double)); + + + /* mark all atoms as unselected */ + + for (i=0; i<m; ++i) { + selected_atoms[i] = 0; + } + + } + + + /* main loop */ + + i=0; + while (resnorm>eps2 && i<T) { + + /* index of next atom */ + memcpy(proj1, proj, m/2*sizeof(double)); + memcpy(proj2, proj + m/2, m/2*sizeof(double)); + for (k=0; k<m/2; k++){ + alpha[k] = (proj1[k] - D1D2[k]*proj2[k])*n12[k]; + beta[k] = (proj2[k] - D1D2[k]*proj1[k])*n12[k]; + } + for (k=0; k<m/2; k++){ + error[k]=0; + for (j=0; j<n; j++){ + error[k]+= (abs(r[j])-D1[k*n+j]*alpha[k]-D2[k*n+j]*beta[k])*(abs(r[j])-D1[k*n+j]*alpha[k]-D2[k*n+j]*beta[k]); + } + } + pos = maxabs(error, m/2); + addproftime(&pd, MAXABS_TIME); + + + /* stop criterion: selected same atom twice, or inner product too small */ + + if (selected_atoms[pos] || alpha[pos]*alpha[pos]<1e-14) { + break; + } + + for (k=0;k<2;k++){ + /* mark selected atom */ + + ind[i] = pos+k*m/2; + selected_atoms[pos+k*m/2] = 1; + + + /* matrix reallocation */ + + if (erroromp && i>=allocated_cols) { + + allocated_cols = (mwSize)(ceil(allocated_cols*MAT_INC_FACTOR) + 1.01); + + Lchol = (double*)mxRealloc(Lchol,n*allocated_cols*sizeof(double)); + + batchomp ? (Gsub = (double*)mxRealloc(Gsub,m*allocated_cols*sizeof(double))) : + (Dsub = (double*)mxRealloc(Dsub,n*allocated_cols*sizeof(double))) ; + } + + + /* append column to Gsub or Dsub */ + + if (batchomp) { + memcpy(Gsub+i*m, G+(pos+k*m/2)*m, m*sizeof(double)); + } + else { + memcpy(Dsub+(i)*n, D+(pos+k*m/2)*n, n*sizeof(double)); + } + + + /*** Cholesky update ***/ + + if (i==0) { + *Lchol = 1; + } + else { + + /* incremental Cholesky decomposition: compute next row of Lchol */ + + if (standardomp) { + matT_vec(1, Dsub, D+n*(pos+k*m/2), tempvec1, n, i); /* compute tempvec1 := Dsub'*d where d is new atom */ + addproftime(&pd, DtD_TIME); + } + else { + vec_assign(tempvec1, Gsub+i*m, ind, i); /* extract tempvec1 := Gsub(ind,i) */ + } + backsubst('L', Lchol, tempvec1, tempvec2, n, i); /* compute tempvec2 = Lchol \ tempvec1 */ + for (j=0; j<i; ++j) { /* write tempvec2 to end of Lchol */ + Lchol[j*n+i] = tempvec2[j]; + } + + /* compute Lchol(i,i) */ + sum = 0; + for (j=0; j<i; ++j) { /* compute sum of squares of last row without Lchol(i,i) */ + sum += SQR(Lchol[j*n+i]); + } + if ( (1-sum) <= 1e-14 ) { /* Lchol(i,i) is zero => selected atoms are dependent */ + break; + } + Lchol[i*n+i] = sqrt(1-sum); + } + + addproftime(&pd, LCHOL_TIME); + + i++; + + } + /* perform orthogonal projection and compute sparse coefficients */ + + vec_assign(tempvec1, DtX + m*signum*DtX_specified, ind, i); /* extract tempvec1 = DtX(ind) */ + cholsolve('L', Lchol, tempvec1, c, n, i); /* solve LL'c = tempvec1 for c */ + addproftime(&pd, COMPCOEF_TIME); + + + /* update alpha = D'*residual */ + + if (standardomp) { + mat_vec(-1, Dsub, c, r, n, i); /* compute r := -Dsub*c */ + vec_sum(1, x+n*signum, r, n); /* compute r := x+r */ + + + /*memcpy(r, x+n*signum, n*sizeof(double)); /* assign r := x */ + /*mat_vec1(-1, Dsub, c, 1, r, n, i); /* compute r := r-Dsub*c */ + + addproftime(&pd, COMPRES_TIME); + matT_vec(1, D, r, proj, n, m); /* compute proj := D'*r */ + addproftime(&pd, DtR_TIME); + + /* update residual norm */ + if (erroromp) { + resnorm = dotprod(r, r, n); + addproftime(&pd, UPDATE_RESNORM_TIME); + } + } + else { + mat_vec(1, Gsub, c, tempvec1, m, i); /* compute tempvec1 := Gsub*c */ + memcpy(proj, DtX + m*signum*DtX_specified, m*sizeof(double)); /* set proj = D'*x */ + vec_sum(-1, tempvec1, proj, m); /* compute proj := proj - tempvec1 */ + addproftime(&pd, UPDATE_DtR_TIME); + + /* update residual norm */ + if (erroromp) { + vec_assign(tempvec2, tempvec1, ind, i); /* assign tempvec2 := tempvec1(ind) */ + delta = dotprod(c,tempvec2,i); /* compute c'*tempvec2 */ + resnorm = resnorm - delta + deltaprev; /* residual norm update */ + deltaprev = delta; + addproftime(&pd, UPDATE_RESNORM_TIME); + } + } + } + + + /*** generate output vector gamma ***/ + + if (gamma_mode == FULL_GAMMA) { /* write the coefs in c to their correct positions in gamma */ + for (j=0; j<i; ++j) { + gammaPr[m*signum + ind[j]] = c[j]; + } + } + else { + /* sort the coefs by index before writing them to gamma */ + quicksort(ind,c,i); + addproftime(&pd, INDEXSORT_TIME); + + /* gamma is full - reallocate */ + if (gamma_count+i >= allocated_coefs) { + + while(gamma_count+i >= allocated_coefs) { + allocated_coefs = (mwSize)(ceil(GAMMA_INC_FACTOR*allocated_coefs) + 1.01); + } + + mxSetNzmax(Gamma, allocated_coefs); + mxSetPr(Gamma, mxRealloc(gammaPr, allocated_coefs*sizeof(double))); + mxSetIr(Gamma, mxRealloc(gammaIr, allocated_coefs*sizeof(mwIndex))); + + gammaPr = mxGetPr(Gamma); + gammaIr = mxGetIr(Gamma); + } + + /* append coefs to gamma and update the indices */ + for (j=0; j<i; ++j) { + gammaPr[gamma_count] = c[j]; + gammaIr[gamma_count] = ind[j]; + gamma_count++; + } + gammaJc[signum+1] = gammaJc[signum] + i; + } + + + + /*** display status messages ***/ + + if (msg_delta>0 && (clock()-lastprint_time)/(double)CLOCKS_PER_SEC >= msg_delta) + { + lastprint_time = clock(); + + /* estimated remainig time */ + secs2hms( ((L-signum-1)/(double)(signum+1)) * ((lastprint_time-starttime)/(double)CLOCKS_PER_SEC) , + &hrs_remain, &mins_remain, &secs_remain); + + mexPrintf("omp: signal %d / %d, estimated remaining time: %02d:%02d:%05.2f\n", + signum+1, L, hrs_remain, mins_remain, secs_remain); + mexEvalString("drawnow;"); + } + + } + + /* end omp */ + + + + /*** print final messages ***/ + + if (msg_delta>0) { + mexPrintf("omp: signal %d / %d\n", signum, L); + } + + if (profile) { + printprofinfo(&pd, erroromp, batchomp, L); + } + + + + /* free memory */ + + if (!DtX_specified) { + mxFree(DtX); + } + if (standardomp) { + mxFree(r); + mxFree(Dsub); + } + else { + mxFree(Gsub); + } + mxFree(tempvec2); + mxFree(tempvec1); + mxFree(Lchol); + mxFree(c); + mxFree(selected_atoms); + mxFree(ind); + mxFree(proj); + mxFree(proj1); + mxFree(proj2); + mxFree(D1); + mxFree(D2); + mxFree(D1D2); + mxFree(n12); + mxFree(alpha); + mxFree(beta); + mxFree(error); + + return Gamma; +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/solvers/SMALL_ompGabor/ompcoreGabor.h Mon Jul 25 17:27:05 2011 +0100 @@ -0,0 +1,80 @@ +/************************************************************************** + * + * File name: ompcore.h + * + * Ron Rubinstein + * Computer Science Department + * Technion, Haifa 32000 Israel + * ronrubin@cs + * + * Last Updated: 18.8.2009 + * + * Contains the core implementation of Batch-OMP / OMP-Cholesky. + * + *************************************************************************/ + + +#ifndef __OMP_CORE_H__ +#define __OMP_CORE_H__ + + +#include "mex.h" + + + +/************************************************************************** + * Perform Batch-OMP or OMP-Cholesky on a specified set of signals, using + * either a fixed number of atoms or an error bound. + * + * Parameters (not all required): + * + * D - the dictionary, of size n X m + * x - the signals, of size n X L + * DtX - D'*x, of size m X L + * XtX - squared norms of the signals in x, sum(x.*x), of length L + * G - D'*D, of size m X m + * T - target sparsity, or maximal number of atoms for error-based OMP + * eps - target residual norm for error-based OMP + * gamma_mode - one of the constants FULL_GAMMA or SPARSE_GAMMA + * profile - if non-zero, profiling info is printed + * msg_delta - positive: the # of seconds between status prints, otherwise: nothing is printed + * erroromp - if nonzero indicates error-based OMP, otherwise fixed sparsity OMP + * + * Usage: + * + * The function can be called using different parameters, and will have + * different complexity depending on the parameters specified. Arrays which + * are not specified should be passed as null (0). When G is specified, + * Batch-OMP is performed. Otherwise, OMP-Cholesky is performed. + * + * Fixed-sparsity usage: + * --------------------- + * Either DtX, or D and x, must be specified. Specifying DtX is more efficient. + * XtX does not need to be specified. + * When D and x are specified, G is not required. However, not providing G + * will significantly degrade efficiency. + * The number of atoms must be specified in T. The value of eps is ignored. + * Finally, set erroromp to 0. + * + * Error-OMP usage: + * ---------------- + * Either DtX and Xtx, or D and x, must be specified. Specifying DtX and XtX + * is more efficient. + * When D and x are specified, G is not required. However, not providing G + * will significantly degrade efficiency. + * The target error must be specified in eps. A hard limit on the number + * of atoms can also be specified via the parameter T. Otherwise, T should + * be negative. Finally, set erroromp to nonzero. + * + * + * Returns: + * An mxArray containing the sparse representations of the signals in x + * (allocated using the appropriate mxCreateXXX() function). + * The array is either full or sparse, depending on gamma_mode. + * + **************************************************************************/ +mxArray* ompcoreGabor(double D[], double x[], double DtX[], double XtX[], double G[], mwSize n, mwSize m, mwSize L, + int T, double eps, int gamma_mode, int profile, double msg_delta, int erroromp); + + +#endif
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/solvers/SMALL_ompGabor/ompmex.c Mon Jul 25 17:27:05 2011 +0100 @@ -0,0 +1,133 @@ +/************************************************************************** + * + * File name: ompmex.c + * + * Ron Rubinstein + * Computer Science Department + * Technion, Haifa 32000 Israel + * ronrubin@cs + * + * Last Updated: 18.8.2009 + * + *************************************************************************/ + +#include "ompcore.h" +#include "omputils.h" +#include "mexutils.h" + + +/* Input Arguments */ + +#define IN_D prhs[0] +#define IN_X prhs[1] +#define IN_DtX prhs[2] +#define IN_G prhs[3] +#define IN_T prhs[4] +#define IN_SPARSE_G prhs[5] +#define IN_MSGDELTA prhs[6] +#define IN_PROFILE prhs[7] + + +/* Output Arguments */ + +#define GAMMA_OUT plhs[0] + + +/***************************************************************************************/ + + +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray*prhs[]) + +{ + double *D, *x, *DtX, *G, msgdelta; + int gmode, profile, T; + mwSize m, n, L; /* D is n x m , X is n x L, DtX is m x L */ + + + /* check parameters */ + + checkmatrix(IN_D, "OMP", "D"); + checkmatrix(IN_X, "OMP", "X"); + checkmatrix(IN_DtX, "OMP", "DtX"); + checkmatrix(IN_G, "OMP", "G"); + + checkscalar(IN_T, "OMP", "T"); + checkscalar(IN_SPARSE_G, "OMP", "sparse_g"); + checkscalar(IN_MSGDELTA, "OMP", "msgdelta"); + checkscalar(IN_PROFILE, "OMP", "profile"); + + + /* get parameters */ + + x = D = DtX = G = 0; + + if (!mxIsEmpty(IN_D)) + D = mxGetPr(IN_D); + + if (!mxIsEmpty(IN_X)) + x = mxGetPr(IN_X); + + if (!mxIsEmpty(IN_DtX)) + DtX = mxGetPr(IN_DtX); + + if (!mxIsEmpty(IN_G)) + G = mxGetPr(IN_G); + + T = (int)(mxGetScalar(IN_T)+1e-2); + if ((int)(mxGetScalar(IN_SPARSE_G)+1e-2)) { + gmode = SPARSE_GAMMA; + } + else { + gmode = FULL_GAMMA; + } + msgdelta = mxGetScalar(IN_MSGDELTA); + profile = (int)(mxGetScalar(IN_PROFILE)+1e-2); + + + /* check sizes */ + + if (D && x) { + n = mxGetM(IN_D); + m = mxGetN(IN_D); + L = mxGetN(IN_X); + + if (mxGetM(IN_X) != n) { + mexErrMsgTxt("D and X have incompatible sizes."); + } + + if (G) { + if (mxGetN(IN_G)!=mxGetM(IN_G)) { + mexErrMsgTxt("G must be a square matrix."); + } + if (mxGetN(IN_G) != m) { + mexErrMsgTxt("D and G have incompatible sizes."); + } + } + } + + else if (DtX) { + m = mxGetM(IN_DtX); + L = mxGetN(IN_DtX); + + n = T; /* arbitrary - it is enough to assume signal length is T */ + + if (mxGetN(IN_G)!=mxGetM(IN_G)) { + mexErrMsgTxt("G must be a square matrix."); + } + if (mxGetN(IN_G) != m) { + mexErrMsgTxt("DtX and G have incompatible sizes."); + } + } + + else { + mexErrMsgTxt("Either D and X, or DtX, must be specified."); + } + + + /* Do OMP! */ + + GAMMA_OUT = ompcore(D, x, DtX, 0, G, n, m, L, T, 0, gmode, profile, msgdelta, 0); + + return; +} +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/solvers/SMALL_ompGabor/ompmex.m Mon Jul 25 17:27:05 2011 +0100 @@ -0,0 +1,22 @@ +%This is the Matlab interface to the OMP MEX implementation. +%The function is not for independent use, only through omp.m. + + +%OMPMEX Matlab interface to the OMP MEX implementation. +% GAMMA = OMPMEX(D,X,DtX,G,L,SPARSE_G,MSGDELTA,PROFILE) invokes the OMP +% MEX function according to the specified parameters. Not all the +% parameters are required. Those among D, X, DtX and G which are not +% specified should be passed as []. +% +% L - the target sparsity. +% SPARSE_G - returns a sparse GAMMA when nonzero, full GAMMA when zero. +% MSGDELTA - the delay in secs between messages. Zero means no messages. +% PROFILE - nonzero means that profiling information should be printed. + + +% Ron Rubinstein +% Computer Science Department +% Technion, Haifa 32000 Israel +% ronrubin@cs +% +% April 2009
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/solvers/SMALL_ompGabor/ompmexGabor.c Mon Jul 25 17:27:05 2011 +0100 @@ -0,0 +1,133 @@ +/************************************************************************** + * + * File name: ompmex.c + * + * Ron Rubinstein + * Computer Science Department + * Technion, Haifa 32000 Israel + * ronrubin@cs + * + * Last Updated: 18.8.2009 + * + *************************************************************************/ + +#include "ompcoreGabor.h" +#include "omputils.h" +#include "mexutils.h" + + +/* Input Arguments */ + +#define IN_D prhs[0] +#define IN_X prhs[1] +#define IN_DtX prhs[2] +#define IN_G prhs[3] +#define IN_T prhs[4] +#define IN_SPARSE_G prhs[5] +#define IN_MSGDELTA prhs[6] +#define IN_PROFILE prhs[7] + + +/* Output Arguments */ + +#define GAMMA_OUT plhs[0] + + +/***************************************************************************************/ + + +void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray*prhs[]) + +{ + double *D, *x, *DtX, *G, msgdelta; + int gmode, profile, T; + mwSize m, n, L; /* D is n x m , X is n x L, DtX is m x L */ + + + /* check parameters */ + + checkmatrix(IN_D, "OMP", "D"); + checkmatrix(IN_X, "OMP", "X"); + checkmatrix(IN_DtX, "OMP", "DtX"); + checkmatrix(IN_G, "OMP", "G"); + + checkscalar(IN_T, "OMP", "T"); + checkscalar(IN_SPARSE_G, "OMP", "sparse_g"); + checkscalar(IN_MSGDELTA, "OMP", "msgdelta"); + checkscalar(IN_PROFILE, "OMP", "profile"); + + + /* get parameters */ + + x = D = DtX = G = 0; + + if (!mxIsEmpty(IN_D)) + D = mxGetPr(IN_D); + + if (!mxIsEmpty(IN_X)) + x = mxGetPr(IN_X); + + if (!mxIsEmpty(IN_DtX)) + DtX = mxGetPr(IN_DtX); + + if (!mxIsEmpty(IN_G)) + G = mxGetPr(IN_G); + + T = (int)(mxGetScalar(IN_T)+1e-2); + if ((int)(mxGetScalar(IN_SPARSE_G)+1e-2)) { + gmode = SPARSE_GAMMA; + } + else { + gmode = FULL_GAMMA; + } + msgdelta = mxGetScalar(IN_MSGDELTA); + profile = (int)(mxGetScalar(IN_PROFILE)+1e-2); + + + /* check sizes */ + + if (D && x) { + n = mxGetM(IN_D); + m = mxGetN(IN_D); + L = mxGetN(IN_X); + + if (mxGetM(IN_X) != n) { + mexErrMsgTxt("D and X have incompatible sizes."); + } + + if (G) { + if (mxGetN(IN_G)!=mxGetM(IN_G)) { + mexErrMsgTxt("G must be a square matrix."); + } + if (mxGetN(IN_G) != m) { + mexErrMsgTxt("D and G have incompatible sizes."); + } + } + } + + else if (DtX) { + m = mxGetM(IN_DtX); + L = mxGetN(IN_DtX); + + n = T; /* arbitrary - it is enough to assume signal length is T */ + + if (mxGetN(IN_G)!=mxGetM(IN_G)) { + mexErrMsgTxt("G must be a square matrix."); + } + if (mxGetN(IN_G) != m) { + mexErrMsgTxt("DtX and G have incompatible sizes."); + } + } + + else { + mexErrMsgTxt("Either D and X, or DtX, must be specified."); + } + + + /* Do OMP! */ + + GAMMA_OUT = ompcoreGabor(D, x, DtX, 0, G, n, m, L, T, 0, gmode, profile, msgdelta, 0); + + return; +} +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/solvers/SMALL_ompGabor/ompmexGabor.m Mon Jul 25 17:27:05 2011 +0100 @@ -0,0 +1,22 @@ +%This is the Matlab interface to the OMP MEX implementation. +%The function is not for independent use, only through omp.m. + + +%OMPMEX Matlab interface to the OMP MEX implementation. +% GAMMA = OMPMEXGABOR(D,X,DtX,G,L,SPARSE_G,MSGDELTA,PROFILE) invokes the OMP +% MEX function according to the specified parameters. Not all the +% parameters are required. Those among D, X, DtX and G which are not +% specified should be passed as []. +% +% L - the target sparsity. +% SPARSE_G - returns a sparse GAMMA when nonzero, full GAMMA when zero. +% MSGDELTA - the delay in secs between messages. Zero means no messages. +% PROFILE - nonzero means that profiling information should be printed. + + +% Ron Rubinstein +% Computer Science Department +% Technion, Haifa 32000 Israel +% ronrubin@cs +% +% April 2009
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/solvers/SMALL_ompGabor/ompprof.c Mon Jul 25 17:27:05 2011 +0100 @@ -0,0 +1,113 @@ +/************************************************************************** + * + * File name: ompprof.c + * + * Ron Rubinstein + * Computer Science Department + * Technion, Haifa 32000 Israel + * ronrubin@cs + * + * Last Updated: 11.4.2009 + * + *************************************************************************/ + + +#include "ompprof.h" + + +/* initialize profiling information */ + +void initprofdata(profdata *pd) +{ + pd->DtX_time = 0; + pd->XtX_time = 0; + pd->DtR_time = 0; + pd->maxabs_time = 0; + pd->DtD_time = 0; + pd->Lchol_time = 0; + pd->compcoef_time = 0; + pd->update_DtR_time = 0; + pd->update_resnorm_time = 0; + pd->compres_time = 0; + pd->indexsort_time = 0; + + pd->DtX_time_counted = 0; + pd->XtX_time_counted = 0; + pd->DtR_time_counted = 0; + pd->DtD_time_counted = 0; + pd->update_DtR_time_counted = 0; + pd->resnorm_time_counted = 0; + pd->compres_time_counted = 0; + pd->indexsort_time_counted = 0; + + pd->prevtime = clock(); +} + + +/* add elapsed time to profiling data according to specified computation */ + +void addproftime(profdata *pd, int comptype) +{ + switch(comptype) { + case DtX_TIME: pd->DtX_time += clock()-pd->prevtime; pd->DtX_time_counted = 1; break; + case XtX_TIME: pd->XtX_time += clock()-pd->prevtime; pd->XtX_time_counted = 1; break; + case DtR_TIME: pd->DtR_time += clock()-pd->prevtime; pd->DtR_time_counted = 1; break; + case DtD_TIME: pd->DtD_time += clock()-pd->prevtime; pd->DtD_time_counted = 1; break; + case COMPRES_TIME: pd->compres_time += clock()-pd->prevtime; pd->compres_time_counted = 1; break; + case UPDATE_DtR_TIME: pd->update_DtR_time += clock()-pd->prevtime; pd->update_DtR_time_counted = 1; break; + case UPDATE_RESNORM_TIME: pd->update_resnorm_time += clock()-pd->prevtime; pd->resnorm_time_counted = 1; break; + case INDEXSORT_TIME: pd->indexsort_time += clock()-pd->prevtime; pd->indexsort_time_counted = 1; break; + case MAXABS_TIME: pd->maxabs_time += clock()-pd->prevtime; break; + case LCHOL_TIME: pd->Lchol_time += clock()-pd->prevtime; break; + case COMPCOEF_TIME: pd->compcoef_time += clock()-pd->prevtime; break; + } + pd->prevtime = clock(); +} + + +/* print profiling info */ + +void printprofinfo(profdata *pd, int erroromp, int batchomp, int signum) +{ + clock_t tottime; + + tottime = pd->DtX_time + pd->XtX_time + pd->DtR_time + pd->DtD_time + pd->compres_time + pd->maxabs_time + + pd->Lchol_time + pd->compcoef_time + pd->update_DtR_time + pd->update_resnorm_time + pd->indexsort_time; + + mexPrintf("\n\n***** Profiling information for %s *****\n\n", erroromp? "OMP2" : "OMP"); + + mexPrintf("OMP mode: %s\n\n", batchomp? "Batch-OMP" : "OMP-Cholesky"); + + mexPrintf("Total signals processed: %d\n\n", signum); + + if (pd->DtX_time_counted) { + mexPrintf("Compute DtX time: %7.3lf seconds\n", pd->DtX_time/(double)CLOCKS_PER_SEC); + } + if (pd->XtX_time_counted) { + mexPrintf("Compute XtX time: %7.3lf seconds\n", pd->XtX_time/(double)CLOCKS_PER_SEC); + } + mexPrintf("Max abs time: %7.3lf seconds\n", pd->maxabs_time/(double)CLOCKS_PER_SEC); + if (pd->DtD_time_counted) { + mexPrintf("Compute DtD time: %7.3lf seconds\n", pd->DtD_time/(double)CLOCKS_PER_SEC); + } + mexPrintf("Lchol update time: %7.3lf seconds\n", pd->Lchol_time/(double)CLOCKS_PER_SEC); + mexPrintf("Compute coef time: %7.3lf seconds\n", pd->compcoef_time/(double)CLOCKS_PER_SEC); + if (pd->compres_time_counted) { + mexPrintf("Compute R time: %7.3lf seconds\n", pd->compres_time/(double)CLOCKS_PER_SEC); + } + if (pd->DtR_time_counted) { + mexPrintf("Compute DtR time: %7.3lf seconds\n", pd->DtR_time/(double)CLOCKS_PER_SEC); + } + if (pd->update_DtR_time_counted) { + mexPrintf("Update DtR time: %7.3lf seconds\n", pd->update_DtR_time/(double)CLOCKS_PER_SEC); + } + if (pd->resnorm_time_counted) { + mexPrintf("Update resnorm time: %7.3lf seconds\n", pd->update_resnorm_time/(double)CLOCKS_PER_SEC); + } + if (pd->indexsort_time_counted) { + mexPrintf("Index sort time: %7.3lf seconds\n", pd->indexsort_time/(double)CLOCKS_PER_SEC); + } + mexPrintf("---------------------------------------\n"); + mexPrintf("Total time: %7.3lf seconds\n\n", tottime/(double)CLOCKS_PER_SEC); +} +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/solvers/SMALL_ompGabor/ompprof.h Mon Jul 25 17:27:05 2011 +0100 @@ -0,0 +1,106 @@ +/************************************************************************** + * + * File name: ompprof.h + * + * Ron Rubinstein + * Computer Science Department + * Technion, Haifa 32000 Israel + * ronrubin@cs + * + * Last Updated: 18.8.2009 + * + * Collection of definitions and functions for profiling the OMP method. + * + *************************************************************************/ + + +#ifndef __OMP_PROF_H__ +#define __OMP_PROF_H__ + +#include "mex.h" +#include <time.h> + + + +/************************************************************************** + * + * Constants and data types. + * + **************************************************************************/ + + +/* constants denoting the various parts of the algorithm */ + +enum { DtX_TIME, XtX_TIME, DtR_TIME, MAXABS_TIME, DtD_TIME, LCHOL_TIME, COMPCOEF_TIME, + UPDATE_DtR_TIME, UPDATE_RESNORM_TIME, COMPRES_TIME, INDEXSORT_TIME }; + + + +/* profiling data container with counters for each part of the algorithm */ + +typedef struct profdata +{ + clock_t prevtime; /* the time when last initialization/call to addproftime() was performed */ + + clock_t DtX_time; + clock_t XtX_time; + clock_t DtR_time; + clock_t maxabs_time; + clock_t DtD_time; + clock_t Lchol_time; + clock_t compcoef_time; + clock_t update_DtR_time; + clock_t update_resnorm_time; + clock_t compres_time; + clock_t indexsort_time; + + /* flags indicating whether profiling data was gathered */ + int DtX_time_counted; + int XtX_time_counted; + int DtR_time_counted; + int DtD_time_counted; + int update_DtR_time_counted; + int resnorm_time_counted; + int compres_time_counted; + int indexsort_time_counted; + +} profdata; + + + +/************************************************************************** + * + * Initialize a profdata structure, zero all counters, and start its timer. + * + **************************************************************************/ +void initprofdata(profdata *pd); + + +/************************************************************************** + * + * Add elapsed time from last call to addproftime(), or from initialization + * of profdata, to the counter specified by comptype. comptype must be one + * of the constants in the enumeration above. + * + **************************************************************************/ +void addproftime(profdata *pd, int comptype); + + +/************************************************************************** + * + * Print the current contents of the counters in profdata. + * + * Parameters: + * pd - the profdata to print + * erroromp - indicates whether error-based (nonzero) or sparsity-based (zero) + * omp was performed. + * batchomp - indicates whether batch-omp (nonzero) or omp-cholesky (zero) + * omp was performed. + * signum - number of signals processed by omp + * + **************************************************************************/ +void printprofinfo(profdata *pd, int erroromp, int batchomp, int signum); + + +#endif +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/solvers/SMALL_ompGabor/omputils.c Mon Jul 25 17:27:05 2011 +0100 @@ -0,0 +1,89 @@ +/************************************************************************** + * + * File name: omputils.c + * + * Ron Rubinstein + * Computer Science Department + * Technion, Haifa 32000 Israel + * ronrubin@cs + * + * Last Updated: 18.8.2009 + * + *************************************************************************/ + +#include "omputils.h" +#include <math.h> + + +const char FULL_GAMMA_STR[] = "full"; +const char SPARSE_GAMMA_STR[] = "sparse"; + + +/* convert seconds to hours, minutes and seconds */ + +void secs2hms(double sectot, int *hrs, int *mins, double *secs) +{ + *hrs = (int)(floor(sectot/3600)+1e-2); + sectot = sectot - 3600*(*hrs); + *mins = (int)(floor(sectot/60)+1e-2); + *secs = sectot - 60*(*mins); +} + + +/* quicksort, public-domain C implementation by Darel Rex Finley. */ +/* modification: sorts the array data[] as well, according to the values in the array vals[] */ + +#define MAX_LEVELS 300 + +void quicksort(mwIndex vals[], double data[], mwIndex n) { + + long piv, beg[MAX_LEVELS], end[MAX_LEVELS], i=0, L, R, swap ; + double datapiv; + + beg[0]=0; + end[0]=n; + + while (i>=0) { + + L=beg[i]; + R=end[i]-1; + + if (L<R) { + + piv=vals[L]; + datapiv=data[L]; + + while (L<R) + { + while (vals[R]>=piv && L<R) + R--; + if (L<R) { + vals[L]=vals[R]; + data[L++]=data[R]; + } + + while (vals[L]<=piv && L<R) + L++; + if (L<R) { + vals[R]=vals[L]; + data[R--]=data[L]; + } + } + + vals[L]=piv; + data[L]=datapiv; + + beg[i+1]=L+1; + end[i+1]=end[i]; + end[i++]=L; + + if (end[i]-beg[i] > end[i-1]-beg[i-1]) { + swap=beg[i]; beg[i]=beg[i-1]; beg[i-1]=swap; + swap=end[i]; end[i]=end[i-1]; end[i-1]=swap; + } + } + else { + i--; + } + } +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/solvers/SMALL_ompGabor/omputils.h Mon Jul 25 17:27:05 2011 +0100 @@ -0,0 +1,77 @@ +/************************************************************************** + * + * File name: omputils.h + * + * Ron Rubinstein + * Computer Science Department + * Technion, Haifa 32000 Israel + * ronrubin@cs + * + * Last Updated: 18.8.2009 + * + * Utility definitions and functions for the OMP library. + * + *************************************************************************/ + + +#ifndef __OMP_UTILS_H__ +#define __OMP_UTILS_H__ + +#include "mex.h" + + +/* constants for the representation mode of gamma */ + +extern const char FULL_GAMMA_STR[]; /* "full" */ +extern const char SPARSE_GAMMA_STR[]; /* "sparse" */ + + +#define FULL_GAMMA 1 +#define SPARSE_GAMMA 2 +#define INVALID_MODE 3 + + + +/************************************************************************** + * Memory management for OMP2. + * + * GAMMA_INC_FACTOR: + * The matrix GAMMA is allocated with sqrt(n)/2 coefficients per signal, + * for a total of nzmax = L*sqrt(n)/2 nonzeros. Whenever GAMMA needs to be + * increased, it is increased by a factor of GAMMA_INC_FACTOR. + * + * MAT_INC_FACTOR: + * The matrices Lchol, Gsub and Dsub are allocated with sqrt(n)/2 + * columns each. If additional columns are needed, this number is + * increased by a factor of MAT_INC_FACTOR. + **************************************************************************/ + +#define GAMMA_INC_FACTOR (1.4) +#define MAT_INC_FACTOR (1.6) + + + +/************************************************************************** + * Convert number of seconds to hour, minute and second representation. + * + * Parameters: + * sectot - total number of seconds + * hrs, mins, secs - output hours (whole) and minutes (whole) and seconds + * + **************************************************************************/ +void secs2hms(double sectot, int *hrs, int *mins, double *secs); + + + +/************************************************************************** + * QuickSort - public-domain C implementation by Darel Rex Finley. + * + * Modified to sort both the array vals[] and the array data[] according + * to the values in the array vals[]. + * + **************************************************************************/ +void quicksort(mwIndex vals[], double data[], mwIndex n); + + +#endif +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/solvers/SMALL_ompGabor/printf.m Mon Jul 25 17:27:05 2011 +0100 @@ -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
--- a/toolboxes/AudioInpaintingToolbox/Experiments/DeclippingExperiment/declipOneSoundExperiment.m Thu Jul 21 16:37:14 2011 +0100 +++ b/toolboxes/AudioInpaintingToolbox/Experiments/DeclippingExperiment/declipOneSoundExperiment.m Mon Jul 25 17:27:05 2011 +0100 @@ -53,7 +53,7 @@ expParam.solver.function = @inpaintSignal_IndependentProcessingOfFrames; expParam.solver.param.N = 512; % frame length expParam.solver.param.N = 256; % frame length - expParam.solver.param.inpaintFrame = @inpaintFrame_OMP_Gabor; % solver function + expParam.solver.param.inpaintFrame = @inpaintFrame_consOMP_Gabor; % solver function expParam.solver.param.OMPerr = 0.001; expParam.solver.param.sparsityDegree = expParam.solver.param.N/4; expParam.solver.param.D_fun = @Gabor_Dictionary; % Dictionary (function handle)
--- a/util/SMALL_init_solver.m Thu Jul 21 16:37:14 2011 +0100 +++ b/util/SMALL_init_solver.m Mon Jul 25 17:27:05 2011 +0100 @@ -1,4 +1,4 @@ -function solver = SMALL_init_solver(varargin) +function solver = SMALL_init_solver(toolbox, name, param, profile) %% Function initialise SMALL structure for sparse representation. % Optional input variables: % toolbox - name of Dictionary Learning toolbox you want to use @@ -17,11 +17,29 @@ % %% -solver.toolbox=[]; -solver.name=[]; -solver.param=[]; -solver.solution=[]; -solver.reconstructed=[]; -solver.time=[]; +if ~ exist( 'toolbox', 'var' ) || isempty(toolbox) + solver.toolbox = []; +else + solver.toolbox = toolbox; +end +if ~ exist( 'name', 'var' ) || isempty(name) + solver.name = []; +else + solver.name = name; +end +if ~ exist( 'param', 'var' ) || isempty(param) + solver.param = []; +else + solver.param = param; +end +if ~ exist( 'profile', 'var' ) || isempty(profile) + solver.profile = 1; +else + solver.profile = profile; +end +solver.add_constraints = 0; +solver.solution = []; +solver.reconstructed = []; +solver.time = []; end \ No newline at end of file
--- a/util/SMALL_solve.m Thu Jul 21 16:37:14 2011 +0100 +++ b/util/SMALL_solve.m Mon Jul 25 17:27:05 2011 +0100 @@ -39,7 +39,10 @@ b = Problem.b; % The right-hand-side vector. end %% -fprintf('\nStarting solver %s... \n', solver.name); +if (solver.profile) + fprintf('\nStarting solver %s... \n', solver.name); +end + start=cputime; tStart=tic; if strcmpi(solver.toolbox,'sparselab') @@ -93,8 +96,10 @@ % Sparse representation time tElapsed=toc(tStart); solver.time = cputime - start; -fprintf('Solver %s finished task in %2f seconds (cpu time). \n', solver.name, solver.time); -fprintf('Solver %s finished task in %2f seconds (tic-toc time). \n', solver.name, tElapsed); +if (solver.profile) + fprintf('Solver %s finished task in %2f seconds (cpu time). \n', solver.name, solver.time); + fprintf('Solver %s finished task in %2f seconds (tic-toc time). \n', solver.name, tElapsed); +end solver.time=tElapsed; % geting around out of memory problem when converting big matrix from % sparse to full... @@ -105,7 +110,7 @@ solver.solution = full(y); end if isfield(Problem,'reconstruct') -% Reconstruct the signal from the solution -solver.reconstructed = Problem.reconstruct(solver.solution); + % Reconstruct the signal from the solution + solver.reconstructed = Problem.reconstruct(solver.solution); end end
--- a/util/ksvd utils/ompbox utils/make.m Thu Jul 21 16:37:14 2011 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,41 +0,0 @@ -function make -%MAKE Build the OMPBox package. -% MAKE compiles all OMPBox MEX functions, using Matlab's default MEX -% compiler. If the MEX compiler has not been set-up before, please run -% -% mex -setup -% -% before using this MAKE file. - -% Ron Rubinstein -% Computer Science Department -% Technion, Haifa 32000 Israel -% ronrubin@cs -% -% August 2009 - - -% detect platform - -compstr = computer; -is64bit = strcmp(compstr(end-1:end),'64'); - - -% compilation parameters - -compile_params = cell(0); -if (is64bit) - compile_params{1} = '-largeArrayDims'; -end - - -% Compile files % - -ompsources = {'mexutils.c','ompcoreGabor.c','omputils.c','myblas.c','ompprof.c'}; - -disp('Compiling ompmex...'); -mex('-g','ompmexGabor.c', ompsources{:},compile_params{:}); - -disp('Compiling omp2mex...'); -mex('-g','omp2mexGabor.c',ompsources{:},compile_params{:}); -
--- a/util/ksvd utils/ompbox utils/mexutils.c Thu Jul 21 16:37:14 2011 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,79 +0,0 @@ -/************************************************************************** - * - * File name: mexutils.c - * - * Ron Rubinstein - * Computer Science Department - * Technion, Haifa 32000 Israel - * ronrubin@cs - * - * Last Updated: 15.8.2009 - * - *************************************************************************/ - -#include "mexutils.h" -#include <math.h> - - - -/* verify that the mxArray contains a double matrix */ - -void checkmatrix(const mxArray *param, char *fname, char *pname) -{ - char errmsg[100]; - sprintf(errmsg, "%.15s requires that %.25s be a double matrix.", fname, pname); - if (!mxIsDouble(param) || mxIsComplex(param) || mxGetNumberOfDimensions(param)>2) { - mexErrMsgTxt(errmsg); - } -} - - -/* verify that the mxArray contains a 1-D double vector */ - -void checkvector(const mxArray *param, char *fname, char *pname) -{ - char errmsg[100]; - sprintf(errmsg, "%.15s requires that %.25s be a double vector.", fname, pname); - if (!mxIsDouble(param) || mxIsComplex(param) || mxGetNumberOfDimensions(param)>2 || (mxGetM(param)!=1 && mxGetN(param)!=1)) { - mexErrMsgTxt(errmsg); - } -} - - -/* verify that the mxArray contains a double scalar */ - -void checkscalar(const mxArray *param, char *fname, char *pname) -{ - char errmsg[100]; - sprintf(errmsg, "%.15s requires that %.25s be a double scalar.", fname, pname); - if (!mxIsDouble(param) || mxIsComplex(param) || mxGetNumberOfDimensions(param)>2 || - mxGetM(param)!=1 || mxGetN(param)!=1) - { - mexErrMsgTxt(errmsg); - } -} - - -/* verify that the mxArray contains a sparse matrix */ - -void checksparse(const mxArray *param, char *fname, char *pname) -{ - char errmsg[100]; - sprintf(errmsg, "%.15s requires that %.25s be sparse.", fname, pname); - if (!mxIsSparse(param)) { - mexErrMsgTxt(errmsg); - } -} - - -/* verify that the mxArray contains a 1-dimensional cell array */ - -void checkcell_1d(const mxArray *param, char *fname, char *pname) -{ - char errmsg[100]; - sprintf(errmsg, "%.15s requires that %.25s be a 1-D cell array.", fname, pname); - if (!mxIsCell(param) || mxGetNumberOfDimensions(param)>2 || (mxGetM(param)!=1 && mxGetN(param)!=1)) { - mexErrMsgTxt(errmsg); - } -} -
--- a/util/ksvd utils/ompbox utils/mexutils.h Thu Jul 21 16:37:14 2011 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,103 +0,0 @@ -/************************************************************************** - * - * File name: mexutils.h - * - * Ron Rubinstein - * Computer Science Department - * Technion, Haifa 32000 Israel - * ronrubin@cs - * - * Last Updated: 18.8.2009 - * - * Utility functions for MEX files. - * - *************************************************************************/ - - -#ifndef __MEX_UTILS_H__ -#define __MEX_UTILS_H__ - -#include "mex.h" - - - -/************************************************************************** - * Function checkmatrix: - * - * Verify that the specified mxArray is real, of type double, and has - * no more than two dimensions. If not, an error message is printed - * and the mex file terminates. - * - * Parameters: - * param - the mxArray to be checked - * fname - the name of the function where the error occured (15 characters or less) - * pname - the name of the parameter (25 characters or less) - * - **************************************************************************/ -void checkmatrix(const mxArray *param, char *fname, char *pname); - - -/************************************************************************** - * Function checkvector: - * - * Verify that the specified mxArray is 1-D, real, and of type double. The - * vector may be a column or row vector. Otherwise, an error message is - * printed and the mex file terminates. - * - * Parameters: - * param - the mxArray to be checked - * fname - the name of the function where the error occured (15 characters or less) - * pname - the name of the parameter (25 characters or less) - * - **************************************************************************/ -void checkvector(const mxArray *param, char *fname, char *pname); - - -/************************************************************************** - * Function checkscalar: - * - * Verify that the specified mxArray represents a real double scalar value. - * If not, an error message is printed and the mex file terminates. - * - * Parameters: - * param - the mxArray to be checked - * fname - the name of the function where the error occured (15 characters or less) - * pname - the name of the parameter (25 characters or less) - * - **************************************************************************/ -void checkscalar(const mxArray *param, char *fname, char *pname); - - -/************************************************************************** - * Function checksparse: - * - * Verify that the specified mxArray contains a sparse matrix. If not, - * an error message is printed and the mex file terminates. - * - * Parameters: - * param - the mxArray to be checked - * fname - the name of the function where the error occured (15 characters or less) - * pname - the name of the parameter (25 characters or less) - * - **************************************************************************/ -void checksparse(const mxArray *param, char *fname, char *pname); - - -/************************************************************************** - * Function checkcell_1d: - * - * Verify that the specified mxArray is a 1-D cell array. The cell array - * may be arranged as either a column or a row. If not, an error message - * is printed and the mex file terminates. - * - * Parameters: - * param - the mxArray to be checked - * fname - the name of the function where the error occured (15 characters or less) - * pname - the name of the parameter (25 characters or less) - * - **************************************************************************/ -void checkcell_1d(const mxArray *param, char *fname, char *pname); - - -#endif -
--- a/util/ksvd utils/ompbox utils/myblas.c Thu Jul 21 16:37:14 2011 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,673 +0,0 @@ -/************************************************************************** - * - * File name: myblas.c - * - * Ron Rubinstein - * Computer Science Department - * Technion, Haifa 32000 Israel - * ronrubin@cs - * - * Version: 1.1 - * Last updated: 13.8.2009 - * - *************************************************************************/ - - -#include "myblas.h" -#include <ctype.h> - - -/* find maximum of absolute values */ - -mwIndex maxabs(double c[], mwSize m) -{ - mwIndex maxid=0, k; - double absval, maxval = SQR(*c); /* use square which is quicker than absolute value */ - - for (k=1; k<m; ++k) { - absval = SQR(c[k]); - if (absval > maxval) { - maxval = absval; - maxid = k; - } - } - return maxid; -} - - -/* compute y := alpha*x + y */ - -void vec_sum(double alpha, double x[], double y[], mwSize n) -{ - mwIndex i; - - for (i=0; i<n; ++i) { - y[i] += alpha*x[i]; - } -} - -/* compute y := alpha*x .* y */ - -void vec_smult(double alpha, double x[], double y[], mwSize n) -{ - mwIndex i; - - for (i=0; i<n; ++i) { - y[i] *= alpha*x[i]; - } -} - -/* compute y := alpha*A*x */ - -void mat_vec(double alpha, double A[], double x[], double y[], mwSize n, mwSize m) -{ - mwIndex i, j, i_n; - double *Ax; - - Ax = mxCalloc(n,sizeof(double)); - - for (i=0; i<m; ++i) { - i_n = i*n; - for (j=0; j<n; ++j) { - Ax[j] += A[i_n+j] * x[i]; - } - } - - for (j=0; j<n; ++j) { - y[j] = alpha*Ax[j]; - } - - mxFree(Ax); -} - - -/* compute y := alpha*A'*x */ - -void matT_vec(double alpha, double A[], double x[], double y[], mwSize n, mwSize m) -{ - mwIndex i, j, n_i; - double sum0, sum1, sum2, sum3; - - for (j=0; j<m; ++j) { - y[j] = 0; - } - - /* use loop unrolling to accelerate computation */ - - for (i=0; i<m; ++i) { - n_i = n*i; - sum0 = sum1 = sum2 = sum3 = 0; - for (j=0; j+4<n; j+=4) { - sum0 += A[n_i+j]*x[j]; - sum1 += A[n_i+j+1]*x[j+1]; - sum2 += A[n_i+j+2]*x[j+2]; - sum3 += A[n_i+j+3]*x[j+3]; - } - y[i] += alpha * ((sum0 + sum1) + (sum2 + sum3)); - while (j<n) { - y[i] += alpha*A[n_i+j]*x[j]; - j++; - } - } -} - - -/* compute y := alpha*A*x */ - -void mat_sp_vec(double alpha, double pr[], mwIndex ir[], mwIndex jc[], double x[], double y[], mwSize n, mwSize m) -{ - - mwIndex i, j, j1, j2; - - for (i=0; i<n; ++i) { - y[i] = 0; - } - - j2 = jc[0]; - for (i=0; i<m; ++i) { - j1 = j2; j2 = jc[i+1]; - for (j=j1; j<j2; ++j) { - y[ir[j]] += alpha * pr[j] * x[i]; - } - } - -} - - -/* compute y := alpha*A'*x */ - -void matT_sp_vec(double alpha, double pr[], mwIndex ir[], mwIndex jc[], double x[], double y[], mwSize n, mwSize m) -{ - - mwIndex i, j, j1, j2; - - for (i=0; i<m; ++i) { - y[i] = 0; - } - - j2 = jc[0]; - for (i=0; i<m; ++i) { - j1 = j2; j2 = jc[i+1]; - for (j=j1; j<j2; ++j) { - y[i] += alpha * pr[j] * x[ir[j]]; - } - } - -} - - -/* compute y := alpha*A*x */ - -void mat_vec_sp(double alpha, double A[], double pr[], mwIndex ir[], mwIndex jc[], double y[], mwSize n, mwSize m) -{ - - mwIndex i, j, j_n, k, kend; - - for (i=0; i<n; ++i) { - y[i] = 0; - } - - kend = jc[1]; - if (kend==0) { /* x is empty */ - return; - } - - for (k=0; k<kend; ++k) { - j = ir[k]; - j_n = j*n; - for (i=0; i<n; ++i) { - y[i] += alpha * A[i+j_n] * pr[k]; - } - } - -} - - -/* compute y := alpha*A'*x */ - -void matT_vec_sp(double alpha, double A[], double pr[], mwIndex ir[], mwIndex jc[], double y[], mwSize n, mwSize m) -{ - - mwIndex i, j, j_n, k, kend; - - for (i=0; i<m; ++i) { - y[i] = 0; - } - - kend = jc[1]; - if (kend==0) { /* x is empty */ - return; - } - - for (j=0; j<m; ++j) { - j_n = j*n; - for (k=0; k<kend; ++k) { - i = ir[k]; - y[j] += alpha * A[i+j_n] * pr[k]; - } - } - -} - - -/* compute y := alpha*A*x */ - -void mat_sp_vec_sp(double alpha, double pr[], mwIndex ir[], mwIndex jc[], double prx[], mwIndex irx[], mwIndex jcx[], double y[], mwSize n, mwSize m) -{ - - mwIndex i, j, k, kend, j1, j2; - - for (i=0; i<n; ++i) { - y[i] = 0; - } - - kend = jcx[1]; - if (kend==0) { /* x is empty */ - return; - } - - for (k=0; k<kend; ++k) { - i = irx[k]; - j1 = jc[i]; j2 = jc[i+1]; - for (j=j1; j<j2; ++j) { - y[ir[j]] += alpha * pr[j] * prx[k]; - } - } - -} - - -/* compute y := alpha*A'*x */ - -void matT_sp_vec_sp(double alpha, double pr[], mwIndex ir[], mwIndex jc[], double prx[], mwIndex irx[], mwIndex jcx[], double y[], mwSize n, mwSize m) -{ - - mwIndex i, j, k, jend, kend, jadd, kadd, delta; - - for (i=0; i<m; ++i) { - y[i] = 0; - } - - kend = jcx[1]; - if (kend==0) { /* x is empty */ - return; - } - - for (i=0; i<m; ++i) { - j = jc[i]; - jend = jc[i+1]; - k = 0; - while (j<jend && k<kend) { - - delta = ir[j] - irx[k]; - - if (delta) { /* if indices differ - increment the smaller one */ - jadd = delta<0; - kadd = 1-jadd; - j += jadd; - k += kadd; - } - - else { /* indices are equal - add to result and increment both */ - y[i] += alpha * pr[j] * prx[k]; - j++; k++; - } - } - } - -} - - -/* matrix-matrix multiplication */ - -void mat_mat(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k) -{ - mwIndex i1, i2, i3, iX, iA, i2_n; - double b; - - for (i1=0; i1<n*k; i1++) { - X[i1] = 0; - } - - for (i2=0; i2<m; ++i2) { - i2_n = i2*n; - iX = 0; - for (i3=0; i3<k; ++i3) { - iA = i2_n; - b = B[i2+i3*m]; - for (i1=0; i1<n; ++i1) { - X[iX++] += A[iA++]*b; - } - } - } - - for (i1=0; i1<n*k; i1++) { - X[i1] *= alpha; - } -} - - -/* matrix-transpose-matrix multiplication */ - -void matT_mat(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k) -{ - mwIndex i1, i2, i3, iX, iA, i2_n; - double *x, sum0, sum1, sum2, sum3; - - for (i2=0; i2<m; ++i2) { - for (i3=0; i3<k; ++i3) { - sum0 = sum1 = sum2 = sum3 = 0; - for (i1=0; i1+4<n; i1+=4) { - sum0 += A[i1+0+i2*n]*B[i1+0+i3*n]; - sum1 += A[i1+1+i2*n]*B[i1+1+i3*n]; - sum2 += A[i1+2+i2*n]*B[i1+2+i3*n]; - sum3 += A[i1+3+i2*n]*B[i1+3+i3*n]; - } - X[i2+i3*m] = (sum0+sum1) + (sum2+sum3); - while(i1<n) { - X[i2+i3*m] += A[i1+i2*n]*B[i1+i3*n]; - i1++; - } - } - } - - for (i1=0; i1<m*k; i1++) { - X[i1] *= alpha; - } -} - - -/* tensor-matrix product */ - -void tens_mat(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k, mwSize l) -{ - mwIndex i1, i2, i3, i4, i2_n, nml; - double b; - - nml = n*m*l; - for (i1=0; i1<nml; ++i1) { - X[i1] = 0; - } - - for (i2=0; i2<m; ++i2) { - i2_n = i2*n; - for (i3=0; i3<k; ++i3) { - for (i4=0; i4<l; ++i4) { - b = B[i4+i3*l]; - for (i1=0; i1<n; ++i1) { - X[i1 + i2_n + i4*n*m] += A[i1 + i2_n + i3*n*m] * b; - } - } - } - } - - for (i1=0; i1<nml; ++i1) { - X[i1] *= alpha; - } -} - - -/* tensor-matrix-transpose product */ - -void tens_matT(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k, mwSize l) -{ - mwIndex i1, i2, i3, i4, i2_n, nml; - double b; - - nml = n*m*l; - for (i1=0; i1<nml; ++i1) { - X[i1] = 0; - } - - for (i2=0; i2<m; ++i2) { - i2_n = i2*n; - for (i4=0; i4<l; ++i4) { - for (i3=0; i3<k; ++i3) { - b = B[i3+i4*k]; - for (i1=0; i1<n; ++i1) { - X[i1 + i2_n + i4*n*m] += A[i1 + i2_n + i3*n*m] * b; - } - } - } - } - - for (i1=0; i1<nml; ++i1) { - X[i1] *= alpha; - } -} - - -/* dot product */ - -double dotprod(double a[], double b[], mwSize n) -{ - double sum = 0; - mwIndex i; - for (i=0; i<n; ++i) - sum += a[i]*b[i]; - return sum; -} - - -/* find maximum of vector */ - -mwIndex maxpos(double c[], mwSize m) -{ - mwIndex maxid=0, k; - double val, maxval = *c; - - for (k=1; k<m; ++k) { - val = c[k]; - if (val > maxval) { - maxval = val; - maxid = k; - } - } - return maxid; -} - - -/* solve L*x = b */ - -void backsubst_L(double L[], double b[], double x[], mwSize n, mwSize k) -{ - mwIndex i, j; - double rhs; - - for (i=0; i<k; ++i) { - rhs = b[i]; - for (j=0; j<i; ++j) { - rhs -= L[j*n+i]*x[j]; - } - x[i] = rhs/L[i*n+i]; - } -} - - -/* solve L'*x = b */ - -void backsubst_Lt(double L[], double b[], double x[], mwSize n, mwSize k) -{ - mwIndex i, j; - double rhs; - - for (i=k; i>=1; --i) { - rhs = b[i-1]; - for (j=i; j<k; ++j) { - rhs -= L[(i-1)*n+j]*x[j]; - } - x[i-1] = rhs/L[(i-1)*n+i-1]; - } -} - - -/* solve U*x = b */ - -void backsubst_U(double U[], double b[], double x[], mwSize n, mwSize k) -{ - mwIndex i, j; - double rhs; - - for (i=k; i>=1; --i) { - rhs = b[i-1]; - for (j=i; j<k; ++j) { - rhs -= U[j*n+i-1]*x[j]; - } - x[i-1] = rhs/U[(i-1)*n+i-1]; - } -} - - -/* solve U'*x = b */ - -void backsubst_Ut(double U[], double b[], double x[], mwSize n, mwSize k) -{ - mwIndex i, j; - double rhs; - - for (i=0; i<k; ++i) { - rhs = b[i]; - for (j=0; j<i; ++j) { - rhs -= U[i*n+j]*x[j]; - } - x[i] = rhs/U[i*n+i]; - } -} - - -/* back substitution solver */ - -void backsubst(char ul, double A[], double b[], double x[], mwSize n, mwSize k) -{ - if (tolower(ul) == 'u') { - backsubst_U(A, b, x, n, k); - } - else if (tolower(ul) == 'l') { - backsubst_L(A, b, x, n, k); - } - else { - mexErrMsgTxt("Invalid triangular matrix type: must be ''U'' or ''L''"); - } -} - - -/* solve equation set using cholesky decomposition */ - -void cholsolve(char ul, double A[], double b[], double x[], mwSize n, mwSize k) -{ - double *tmp; - - tmp = mxMalloc(k*sizeof(double)); - - if (tolower(ul) == 'l') { - backsubst_L(A, b, tmp, n, k); - backsubst_Lt(A, tmp, x, n, k); - } - else if (tolower(ul) == 'u') { - backsubst_Ut(A, b, tmp, n, k); - backsubst_U(A, tmp, x, n, k); - } - else { - mexErrMsgTxt("Invalid triangular matrix type: must be either ''U'' or ''L''"); - } - - mxFree(tmp); -} - - -/* perform a permutation assignment y := x(ind(1:k)) */ - -void vec_assign(double y[], double x[], mwIndex ind[], mwSize k) -{ - mwIndex i; - - for (i=0; i<k; ++i) - y[i] = x[ind[i]]; -} - - -/* matrix transpose */ - -void transpose(double X[], double Y[], mwSize n, mwSize m) -{ - mwIndex i, j, i_m, j_n; - - if (n<m) { - for (j=0; j<m; ++j) { - j_n = j*n; - for (i=0; i<n; ++i) { - Y[j+i*m] = X[i+j_n]; - } - } - } - else { - for (i=0; i<n; ++i) { - i_m = i*m; - for (j=0; j<m; ++j) { - Y[j+i_m] = X[i+j*n]; - } - } - } -} - - -/* print contents of matrix */ - -void printmat(double A[], int n, int m, char* matname) -{ - int i, j; - mexPrintf("\n%s = \n\n", matname); - - if (n*m==0) { - mexPrintf(" Empty matrix: %d-by-%d\n\n", n, m); - return; - } - - for (i=0; i<n; ++i) { - for (j=0; j<m; ++j) - mexPrintf(" %lf", A[j*n+i]); - mexPrintf("\n"); - } - mexPrintf("\n"); -} - - -/* print contents of sparse matrix */ - -void printspmat(mxArray *a, char* matname) -{ - mwIndex *aJc = mxGetJc(a); - mwIndex *aIr = mxGetIr(a); - double *aPr = mxGetPr(a); - - int i; - - mexPrintf("\n%s = \n\n", matname); - - for (i=0; i<aJc[1]; ++i) - printf(" (%d,1) = %lf\n", aIr[i]+1,aPr[i]); - - mexPrintf("\n"); -} - - - -/* matrix multiplication using Winograd's algorithm */ - -/* -void mat_mat2(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k) -{ - - mwIndex i1, i2, i3, iX, iA, i2_n; - double b, *AA, *BB; - - AA = mxCalloc(n,sizeof(double)); - BB = mxCalloc(k,sizeof(double)); - - for (i1=0; i1<n*k; i1++) { - X[i1] = 0; - } - - for (i1=0; i1<n; ++i1) { - for (i2=0; i2<m/2; ++i2) { - AA[i1] += A[i1+2*i2*n]*A[i1+(2*i2+1)*n]; - } - } - - for (i2=0; i2<k; ++i2) { - for (i1=0; i1<m/2; ++i1) { - BB[i2] += B[2*i1+i2*m]*B[2*i1+1+i2*m]; - } - } - - for (i2=0; i2<k; ++i2) { - for (i3=0; i3<m/2; ++i3) { - for (i1=0; i1<n; ++i1) { - X[i1+i2*n] += (A[i1+(2*i3)*n]+B[2*i3+1+i2*m])*(A[i1+(2*i3+1)*n]+B[2*i3+i2*m]); - } - } - } - - if (m%2) { - for (i2=0; i2<k; ++i2) { - for (i1=0; i1<n; ++i1) { - X[i1+i2*n] += A[i1+(m-1)*n]*B[m-1+i2*m]; - } - } - } - - for (i2=0; i2<k; ++i2) { - for (i1=0; i1<n; ++i1) { - X[i1+i2*n] -= (AA[i1] + BB[i2]); - X[i1+i2*n] *= alpha; - } - } - - mxFree(AA); - mxFree(BB); -} -*/ - - - -
--- a/util/ksvd utils/ompbox utils/myblas.h Thu Jul 21 16:37:14 2011 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,511 +0,0 @@ -/************************************************************************** - * - * File name: myblas.h - * - * Ron Rubinstein - * Computer Science Department - * Technion, Haifa 32000 Israel - * ronrubin@cs - * - * Version: 1.1 - * Last updated: 17.8.2009 - * - * A collection of basic linear algebra functions, in the spirit of the - * BLAS/LAPACK libraries. - * - *************************************************************************/ - - - -#ifndef __MY_BLAS_H__ -#define __MY_BLAS_H__ - - -#include "mex.h" -#include <math.h> - - - -/************************************************************************** - * Squared value. - **************************************************************************/ -#define SQR(X) ((X)*(X)) - - - -/************************************************************************** - * Matrix-vector multiplication. - * - * Computes an operation of the form: - * - * y := alpha*A*x - * - * Parameters: - * A - matrix of size n X m - * x - vector of length m - * y - output vector of length n - * alpha - real constant - * n, m - dimensions of A - * - * Note: This function re-writes the contents of y. - * - **************************************************************************/ -void mat_vec(double alpha, double A[], double x[], double y[], mwSize n, mwSize m); - - - -/************************************************************************** - * Matrix-transpose-vector multiplication. - * - * Computes an operation of the form: - * - * y := alpha*A'*x - * - * Parameters: - * A - matrix of size n X m - * x - vector of length n - * y - output vector of length m - * alpha - real constant - * n, m - dimensions of A - * - * Note: This function re-writes the contents of y. - * - **************************************************************************/ -void matT_vec(double alpha, double A[], double x[], double y[], mwSize n, mwSize m); - - - -/************************************************************************** - * Sparse-matrix-vector multiplication. - * - * Computes an operation of the form: - * - * y := alpha*A*x - * - * where A is a sparse matrix. - * - * Parameters: - * pr,ir,jc - sparse representation of the matrix A, of size n x m - * x - vector of length m - * y - output vector of length n - * alpha - real constant - * n, m - dimensions of A - * - * Note: This function re-writes the contents of y. - * - **************************************************************************/ -void mat_sp_vec(double alpha, double pr[], mwIndex ir[], mwIndex jc[], double x[], double y[], mwSize n, mwSize m); - - - -/************************************************************************** - * Sparse-matrix-transpose-vector multiplication. - * - * Computes an operation of the form: - * - * y := alpha*A'*x - * - * where A is a sparse matrix. - * - * Parameters: - * pr,ir,jc - sparse representation of the matrix A, of size n x m - * x - vector of length m - * y - output vector of length n - * alpha - real constant - * n, m - dimensions of A - * - * Note: This function re-writes the contents of y. - * - **************************************************************************/ -void matT_sp_vec(double alpha, double pr[], mwIndex ir[], mwIndex jc[], double x[], double y[], mwSize n, mwSize m); - - - -/************************************************************************** - * Matrix-sparse-vector multiplication. - * - * Computes an operation of the form: - * - * y := alpha*A*x - * - * where A is a matrix and x is a sparse vector. - * - * Parameters: - * A - matrix of size n X m - * pr,ir,jc - sparse representation of the vector x, of length m - * y - output vector of length n - * alpha - real constant - * n, m - dimensions of A - * - * Note: This function re-writes the contents of y. - * - **************************************************************************/ -void mat_vec_sp(double alpha, double A[], double pr[], mwIndex ir[], mwIndex jc[], double y[], mwSize n, mwSize m); - - - -/************************************************************************** - * Matrix-transpose-sparse-vector multiplication. - * - * Computes an operation of the form: - * - * y := alpha*A'*x - * - * where A is a matrix and x is a sparse vector. - * - * Parameters: - * A - matrix of size n X m - * pr,ir,jc - sparse representation of the vector x, of length n - * y - output vector of length m - * alpha - real constant - * n, m - dimensions of A - * - * Note: This function re-writes the contents of y. - * - **************************************************************************/ -void matT_vec_sp(double alpha, double A[], double pr[], mwIndex ir[], mwIndex jc[], double y[], mwSize n, mwSize m); - - - -/************************************************************************** - * Sparse-matrix-sparse-vector multiplication. - * - * Computes an operation of the form: - * - * y := alpha*A*x - * - * where A is a sparse matrix and x is a sparse vector. - * - * Parameters: - * pr,ir,jc - sparse representation of the matrix A, of size n x m - * prx,irx,jcx - sparse representation of the vector x (of length m) - * y - output vector of length n - * alpha - real constant - * n, m - dimensions of A - * - * Note: This function re-writes the contents of y. - * - **************************************************************************/ -void mat_sp_vec_sp(double alpha, double pr[], mwIndex ir[], mwIndex jc[], double prx[], mwIndex irx[], mwIndex jcx[], double y[], mwSize n, mwSize m); - - - -/************************************************************************** - * Sparse-matrix-transpose-sparse-vector multiplication. - * - * Computes an operation of the form: - * - * y := alpha*A'*x - * - * where A is a sparse matrix and x is a sparse vector. - * - * Importnant note: this function is provided for completeness, but is NOT efficient. - * If possible, convert x to non-sparse representation and use matT_vec_sp instead. - * - * Parameters: - * pr,ir,jc - sparse representation of the matrix A, of size n x m - * prx,irx,jcx - sparse representation of the vector x (of length n) - * y - output vector of length n - * alpha - real constant - * n, m - dimensions of A - * - * Note: This function re-writes the contents of y. - * - **************************************************************************/ -void matT_sp_vec_sp(double alpha, double pr[], mwIndex ir[], mwIndex jc[], double prx[], mwIndex irx[], mwIndex jcx[], double y[], mwSize n, mwSize m); - - - -/************************************************************************** - * Matrix-matrix multiplication. - * - * Computes an operation of the form: - * - * X := alpha*A*B - * - * Parameters: - * A - matrix of size n X m - * B - matrix of size m X k - * X - output matrix of size n X k - * alpha - real constant - * n, m, k - dimensions of A, B - * - * Note: This function re-writes the contents of X. - * - **************************************************************************/ -void mat_mat(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k); - - - -/************************************************************************** - * Matrix-transpose-matrix multiplication. - * - * Computes an operation of the form: - * - * X := alpha*A*B - * - * Parameters: - * A - matrix of size n X m - * B - matrix of size m X k - * X - output matrix of size n X k - * alpha - real constant - * n, m, k - dimensions of A, B - * - * Note: This function re-writes the contents of X. - * - **************************************************************************/ -void matT_mat(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k); - - - -/************************************************************************** - * Tensor-matrix multiplication. - * - * This function accepts a 3-D tensor A of size n X m X k - * and a 2-D matrix B of size l X k. - * The function computes the 3-D tensor X of size n X m X l, where - * - * X(i,j,:) = B*A(i,j,:) - * - * for all i,j. - * - * Parameters: - * A - tensor of size n X m X k - * B - matrix of size l X k - * X - output tensor of size n X m X l - * alpha - real constant - * n, m, k, l - dimensions of A, B - * - * Note: This function re-writes the contents of X. - * - **************************************************************************/ -void tens_mat(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k, mwSize l); - - - -/************************************************************************** - * Tensor-matrix-transpose multiplication. - * - * This function accepts a 3-D tensor A of size n X m X k - * and a 2-D matrix B of size k X l. - * The function computes the 3-D tensor X of size n X m X l, where - * - * X(i,j,:) = B'*A(i,j,:) - * - * for all i,j. - * - * Parameters: - * A - tensor of size n X m X k - * B - matrix of size k X l - * X - output tensor of size n X m X l - * alpha - real constant - * n, m, k, l - dimensions of A, B - * - * Note: This function re-writes the contents of X. - * - **************************************************************************/ -void tens_matT(double alpha, double A[], double B[], double X[], mwSize n, mwSize m, mwSize k, mwSize l); - - - -/************************************************************************** - * Vector-vector sum. - * - * Computes an operation of the form: - * - * y := alpha*x + y - * - * Parameters: - * x - vector of length n - * y - output vector of length n - * alpha - real constant - * n - length of x,y - * - * Note: This function re-writes the contents of y. - * - **************************************************************************/ -void vec_sum(double alpha, double x[], double y[], mwSize n); - -/************************************************************************** - * Vector-vector scalar multiply. - * - * Computes an operation of the form: - * - * y := alpha* x.*y - * - * Parameters: - * x - vector of length n - * y - output vector of length n - * alpha - real constant - * n - length of x,y - * - * Note: This function re-writes the contents of y. - * - **************************************************************************/ - - -void vec_smult(double alpha, double x[], double y[], mwSize n); - - -/************************************************************************** - * Triangular back substitution. - * - * Solve the set of linear equations - * - * T*x = b - * - * where T is lower or upper triangular. - * - * Parameters: - * ul - 'U' for upper triangular, 'L' for lower triangular - * A - matrix of size n x m containing T - * b - vector of length k - * x - output vector of length k - * n - size of first dimension of A - * k - the size of the equation set, k<=n,m - * - * Note: - * The matrix A can be of any size n X m, as long as n,m >= k. - * Only the lower/upper triangle of the submatrix A(1:k,1:k) defines the - * matrix T (depending on the parameter ul). - * - **************************************************************************/ -void backsubst(char ul, double A[], double b[], double x[], mwSize n, mwSize k); - - - -/************************************************************************** - * Solve a set of equations using a Cholesky decomposition. - * - * Solve the set of linear equations - * - * M*x = b - * - * where M is positive definite with a known Cholesky decomposition: - * either M=L*L' (L lower triangular) or M=U'*U (U upper triangular). - * - * Parameters: - * ul - 'U' for upper triangular, 'L' for lower triangular decomposition - * A - matrix of size n x m with the Cholesky decomposition of M - * b - vector of length k - * x - output vector of length k - * n - size of first dimension of A - * k - the size of the equation set, k<=n,m - * - * Note: - * The matrix A can be of any size n X m, as long as n,m >= k. - * Only the lower/upper triangle of the submatrix A(1:k,1:k) is used as - * the Cholesky decomposition of M (depending on the parameter ul). - * - **************************************************************************/ -void cholsolve(char ul, double A[], double b[], double x[], mwSize n, mwSize k); - - - -/************************************************************************** - * Maximum absolute value. - * - * Returns the index of the coefficient with maximal absolute value in a vector. - * - * Parameters: - * x - vector of length n - * n - length of x - * - **************************************************************************/ -mwIndex maxabs(double x[], mwSize n); - - - -/************************************************************************** - * Maximum vector element. - * - * Returns the index of the maximal coefficient in a vector. - * - * Parameters: - * x - vector of length n - * n - length of x - * - **************************************************************************/ -mwIndex maxpos(double x[], mwSize n); - - - -/************************************************************************** - * Vector-vector dot product. - * - * Computes an operation of the form: - * - * c = a'*b - * - * Parameters: - * a, b - vectors of length n - * n - length of a,b - * - * Returns: The dot product c. - * - **************************************************************************/ -double dotprod(double a[], double b[], mwSize n); - - - -/************************************************************************** - * Indexed vector assignment. - * - * Perform a permutation assignment of the form - * - * y = x(ind) - * - * where ind is an array of indices to x. - * - * Parameters: - * y - output vector of length k - * x - input vector of arbitrary length - * ind - array of indices into x (indices begin at 0) - * k - length of the array ind - * - **************************************************************************/ -void vec_assign(double y[], double x[], mwIndex ind[], mwSize k); - - - -/************************************************************************** - * Matrix transpose. - * - * Computes Y := X' - * - * Parameters: - * X - input matrix of size n X m - * Y - output matrix of size m X n - * n, m - dimensions of X - * - **************************************************************************/ -void transpose(double X[], double Y[], mwSize n, mwSize m); - - - -/************************************************************************** - * Print a matrix. - * - * Parameters: - * A - matrix of size n X m - * n, m - dimensions of A - * matname - name of matrix to display - * - **************************************************************************/ -void printmat(double A[], int n, int m, char* matname); - - - -/************************************************************************** - * Print a sparse matrix. - * - * Parameters: - * A - sparse matrix of type double - * matname - name of matrix to display - * - **************************************************************************/ -void printspmat(mxArray *A, char* matname); - - -#endif -
--- a/util/ksvd utils/ompbox utils/omp2Gabor.m Thu Jul 21 16:37:14 2011 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,206 +0,0 @@ -function gamma = omp2(varargin) -%OMP2 Error-constrained Orthogonal Matching Pursuit. -% GAMMA = OMP2(D,X,G,EPSILON) solves the optimization problem -% -% min |GAMMA|_0 s.t. |X - D*GAMMA|_2 <= EPSILON -% gamma -% -% for each of the signals in X, using Batch Orthogonal Matching Pursuit. -% Here, D is a dictionary with normalized columns, X is a matrix -% containing column signals, EPSILON is the error target for each signal, -% and G is the Gramm matrix D'*D. The output GAMMA is a matrix containing -% the sparse representations as its columns. -% -% GAMMA = OMP2(D,X,[],EPSILON) performs the same operation, but without -% the matrix G, using OMP-Cholesky. This call produces the same output as -% Batch-OMP, but is significantly slower. Using this syntax is only -% recommended when available memory is too small to store G. -% -% GAMMA = OMP2(DtX,XtX,G,EPSILON) is the fastest implementation of OMP2, -% but also requires the most memory. Here, DtX stores the projections -% D'*X, and XtX is a row vector containing the squared norms of the -% signals, sum(X.*X). In this case Batch-OMP is used, but without having -% to compute D'*X and XtX in advance, which slightly improves runtime. -% Note that in general, the call -% -% GAMMA = OMP2(D'*X, sum(X.*X), G, EPSILON); -% -% will be faster than the call -% -% GAMMA = OMP2(D,X,G,EPSILON); -% -% due to optimized matrix multiplications in Matlab. However, when the -% entire matrix D'*X cannot be stored in memory, one of the other two -% versions can be used. Both compute D'*X for just one signal at a time, -% and thus require much less memory. -% -% GAMMA = OMP2(...,PARAM1,VAL1,PARAM2,VAL2,...) specifies additional -% parameters for OMP2. Available parameters are: -% -% 'gammamode' - Specifies the representation mode for GAMMA. Can be -% either 'full' or 'sparse', corresponding to a full or -% sparse matrix, respectively. By default, GAMMA is -% returned as a sparse matrix. -% 'maxatoms' - Limits the number of atoms in the representation of each -% signal. If specified, the number of atoms in each -% representation does not exceed this number, even if the -% error target is not met. Specifying maxatoms<0 implies -% no limit (default). -% 'messages' - Specifies whether progress messages should be displayed. -% When positive, this is the number of seconds between -% status prints. When negative, indicates that no messages -% should be displayed (this is the default). -% 'checkdict' - Specifies whether dictionary normalization should be -% verified. When set to 'on' (default) the dictionary -% atoms are verified to be of unit L2-norm. Setting this -% parameter to 'off' disables verification and accelerates -% function performance. Note that an unnormalized -% dictionary will produce invalid results. -% 'profile' - Can be either 'on' or 'off'. When 'on', profiling -% information is displayed at the end of the funciton -% execution. -% -% -% Summary of OMP2 versions: -% -% version | speed | memory -% ------------------------------------------------------------- -% OMP2(DtX,XtX,G,EPSILON) | very fast | very large -% OMP2(D,X,G,EPSILON) | fast | moderate -% OMP2(D,X,[],EPSILON) | very slow | small -% ------------------------------------------------------------- -% -% -% References: -% [1] M. Elad, R. Rubinstein, and M. Zibulevsky, "Efficient Implementation -% of the K-SVD Algorithm using Batch Orthogonal Matching Pursuit", -% Technical Report - CS, Technion, April 2008. -% -% See also OMP. - - -% Ron Rubinstein -% Computer Science Department -% Technion, Haifa 32000 Israel -% ronrubin@cs -% -% April 2009 - - -% default options - -sparse_gamma = 1; -msgdelta = -1; -maxatoms = -1; -checkdict = 1; -profile = 0; - - -% determine number of parameters - -paramnum = 1; -while (paramnum<=nargin && ~ischar(varargin{paramnum})) - paramnum = paramnum+1; -end -paramnum = paramnum-1; - - -% parse options - -for i = paramnum+1:2:length(varargin) - paramname = varargin{i}; - paramval = varargin{i+1}; - - switch lower(paramname) - - case 'gammamode' - if (strcmpi(paramval,'sparse')) - sparse_gamma = 1; - elseif (strcmpi(paramval,'full')) - sparse_gamma = 0; - else - error('Invalid GAMMA mode'); - end - - case 'maxatoms' - maxatoms = paramval; - - case 'messages' - msgdelta = paramval; - - case 'checkdict' - if (strcmpi(paramval,'on')) - checkdict = 1; - elseif (strcmpi(paramval,'off')) - checkdict = 0; - else - error('Invalid checkdict option'); - end - - case 'profile' - if (strcmpi(paramval,'on')) - profile = 1; - elseif (strcmpi(paramval,'off')) - profile = 0; - else - error('Invalid profile mode'); - end - - otherwise - error(['Unknown option: ' paramname]); - end - -end - - -% determine call type - -if (paramnum==4) - - n1 = size(varargin{1},1); - n2 = size(varargin{2},1); - n3 = size(varargin{3},1); - - if ( (n1>1 && n2==1) || (n1==1 && n2==1 && n3==1) ) % DtX,XtX,G,EPSILON - - DtX = varargin{1}; - XtX = varargin{2}; - G = varargin{3}; - epsilon = varargin{4}; - D = []; - X = []; - - else % D,X,G,EPSILON - - D = varargin{1}; - X = varargin{2}; - G = varargin{3}; - epsilon = varargin{4}; - DtX = []; - XtX = []; - - end - -else - error('Invalid number of parameters'); -end - -G=[]; - -% verify dictionary normalization - -if (checkdict) - if (isempty(G)) - atomnorms = sum(D.*D); - else - atomnorms = diag(G); - end - if (any(abs(atomnorms-1) > 1e-2)) - error('Dictionary columns must be normalized to unit length'); - end -end - - -% omp - -gamma = omp2mexGabor(D,X,DtX,XtX,G,epsilon,sparse_gamma,msgdelta,maxatoms,profile);
--- a/util/ksvd utils/ompbox utils/omp2mex.c Thu Jul 21 16:37:14 2011 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,156 +0,0 @@ -/************************************************************************** - * - * File name: omp2mex.c - * - * Ron Rubinstein - * Computer Science Department - * Technion, Haifa 32000 Israel - * ronrubin@cs - * - * Last Updated: 18.8.2009 - * - *************************************************************************/ - -#include "ompcore.h" -#include "omputils.h" -#include "mexutils.h" - - -/* Input Arguments */ - -#define IN_D prhs[0] -#define IN_X prhs[1] -#define IN_DtX prhs[2] -#define IN_XtX prhs[3] -#define IN_G prhs[4] -#define IN_EPS prhs[5] -#define IN_SPARSE_G prhs[6] -#define IN_MSGDELTA prhs[7] -#define IN_MAXATOMS prhs[8] -#define IN_PROFILE prhs[9] - - -/* Output Arguments */ - -#define GAMMA_OUT plhs[0] - - -/***************************************************************************************/ - - -void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray*prhs[]) - -{ - double *D, *x, *DtX, *XtX, *G, eps, msgdelta; - int gmode, maxatoms, profile; - mwSize m, n, L; /* D is n x m , X is n x L, DtX is m x L */ - - - /* check parameters */ - - checkmatrix(IN_D, "OMP2", "D"); - checkmatrix(IN_X, "OMP2", "X"); - checkmatrix(IN_DtX, "OMP2", "DtX"); - checkmatrix(IN_XtX, "OMP2", "XtX"); - checkmatrix(IN_G, "OMP2", "G"); - - checkscalar(IN_EPS, "OMP2", "EPSILON"); - checkscalar(IN_SPARSE_G, "OMP2", "sparse_g"); - checkscalar(IN_MSGDELTA, "OMP2", "msgdelta"); - checkscalar(IN_MAXATOMS, "OMP2", "maxatoms"); - checkscalar(IN_PROFILE, "OMP2", "profile"); - - - /* get parameters */ - - x = D = DtX = XtX = G = 0; - - if (!mxIsEmpty(IN_D)) - D = mxGetPr(IN_D); - - if (!mxIsEmpty(IN_X)) - x = mxGetPr(IN_X); - - if (!mxIsEmpty(IN_DtX)) - DtX = mxGetPr(IN_DtX); - - if (!mxIsEmpty(IN_XtX)) - XtX = mxGetPr(IN_XtX); - - if (!mxIsEmpty(IN_G)) - G = mxGetPr(IN_G); - - eps = mxGetScalar(IN_EPS); - if ((int)(mxGetScalar(IN_SPARSE_G)+1e-2)) { - gmode = SPARSE_GAMMA; - } - else { - gmode = FULL_GAMMA; - } - msgdelta = mxGetScalar(IN_MSGDELTA); - if (mxGetScalar(IN_MAXATOMS) < -1e-5) { - maxatoms = -1; - } - else { - maxatoms = (int)(mxGetScalar(IN_MAXATOMS)+1e-2); - } - profile = (int)(mxGetScalar(IN_PROFILE)+1e-2); - - - /* check sizes */ - - if (D && x) { - n = mxGetM(IN_D); - m = mxGetN(IN_D); - L = mxGetN(IN_X); - - if (mxGetM(IN_X) != n) { - mexErrMsgTxt("D and X have incompatible sizes."); - } - - if (G) { - if (mxGetN(IN_G)!=mxGetM(IN_G)) { - mexErrMsgTxt("G must be a square matrix."); - } - if (mxGetN(IN_G) != m) { - mexErrMsgTxt("D and G have incompatible sizes."); - } - } - } - - else if (DtX && XtX) { - m = mxGetM(IN_DtX); - L = mxGetN(IN_DtX); - - /* set n to an arbitrary value that is at least the max possible number of selected atoms */ - - if (maxatoms>0) { - n = maxatoms; - } - else { - n = m; - } - - if ( !(mxGetM(IN_XtX)==L && mxGetN(IN_XtX)==1) && !(mxGetM(IN_XtX)==1 && mxGetN(IN_XtX)==L) ) { - mexErrMsgTxt("DtX and XtX have incompatible sizes."); - } - - if (mxGetN(IN_G)!=mxGetM(IN_G)) { - mexErrMsgTxt("G must be a square matrix."); - } - if (mxGetN(IN_G) != m) { - mexErrMsgTxt("DtX and G have incompatible sizes."); - } - } - - else { - mexErrMsgTxt("Either D and X, or DtX and XtX, must be specified."); - } - - - /* Do OMP! */ - - GAMMA_OUT = ompcore(D, x, DtX, XtX, G, n, m, L, maxatoms, eps, gmode, profile, msgdelta, 1); - - return; -}
--- a/util/ksvd utils/ompbox utils/omp2mex.m Thu Jul 21 16:37:14 2011 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,23 +0,0 @@ -%This is the Matlab interface to the OMP2 MEX implementation. -%The function is not for independent use, only through omp2.m. - - -%OMP2MEX Matlab interface to the OMP2 MEX implementation. -% GAMMA = OMP2MEX(D,X,DtX,XtX,G,EPSILON,SPARSE_G,MSGDELTA,MAXATOMS,PROFILE) -% invokes the OMP2 MEX function according to the specified parameters. Not -% all the parameters are required. Those among D, X, DtX, XtX and G which -% are not specified should be passed as []. -% -% EPSILON - the target error. -% SPARSE_G - returns a sparse GAMMA when nonzero, full GAMMA when zero. -% MSGDELTA - the delay in secs between messages. Zero means no messages. -% MAXATOMS - the max number of atoms per signal, negative for no max. -% PROFILE - nonzero means that profiling information should be printed. - - -% Ron Rubinstein -% Computer Science Department -% Technion, Haifa 32000 Israel -% ronrubin@cs -% -% April 2009
--- a/util/ksvd utils/ompbox utils/omp2mexGabor.c Thu Jul 21 16:37:14 2011 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,156 +0,0 @@ -/************************************************************************** - * - * File name: omp2mex.c - * - * Ron Rubinstein - * Computer Science Department - * Technion, Haifa 32000 Israel - * ronrubin@cs - * - * Last Updated: 18.8.2009 - * - *************************************************************************/ - -#include "ompcoreGabor.h" -#include "omputils.h" -#include "mexutils.h" - - -/* Input Arguments */ - -#define IN_D prhs[0] -#define IN_X prhs[1] -#define IN_DtX prhs[2] -#define IN_XtX prhs[3] -#define IN_G prhs[4] -#define IN_EPS prhs[5] -#define IN_SPARSE_G prhs[6] -#define IN_MSGDELTA prhs[7] -#define IN_MAXATOMS prhs[8] -#define IN_PROFILE prhs[9] - - -/* Output Arguments */ - -#define GAMMA_OUT plhs[0] - - -/***************************************************************************************/ - - -void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray*prhs[]) - -{ - double *D, *x, *DtX, *XtX, *G, eps, msgdelta; - int gmode, maxatoms, profile; - mwSize m, n, L; /* D is n x m , X is n x L, DtX is m x L */ - - - /* check parameters */ - - checkmatrix(IN_D, "OMP2", "D"); - checkmatrix(IN_X, "OMP2", "X"); - checkmatrix(IN_DtX, "OMP2", "DtX"); - checkmatrix(IN_XtX, "OMP2", "XtX"); - checkmatrix(IN_G, "OMP2", "G"); - - checkscalar(IN_EPS, "OMP2", "EPSILON"); - checkscalar(IN_SPARSE_G, "OMP2", "sparse_g"); - checkscalar(IN_MSGDELTA, "OMP2", "msgdelta"); - checkscalar(IN_MAXATOMS, "OMP2", "maxatoms"); - checkscalar(IN_PROFILE, "OMP2", "profile"); - - - /* get parameters */ - - x = D = DtX = XtX = G = 0; - - if (!mxIsEmpty(IN_D)) - D = mxGetPr(IN_D); - - if (!mxIsEmpty(IN_X)) - x = mxGetPr(IN_X); - - if (!mxIsEmpty(IN_DtX)) - DtX = mxGetPr(IN_DtX); - - if (!mxIsEmpty(IN_XtX)) - XtX = mxGetPr(IN_XtX); - - if (!mxIsEmpty(IN_G)) - G = mxGetPr(IN_G); - - eps = mxGetScalar(IN_EPS); - if ((int)(mxGetScalar(IN_SPARSE_G)+1e-2)) { - gmode = SPARSE_GAMMA; - } - else { - gmode = FULL_GAMMA; - } - msgdelta = mxGetScalar(IN_MSGDELTA); - if (mxGetScalar(IN_MAXATOMS) < -1e-5) { - maxatoms = -1; - } - else { - maxatoms = (int)(mxGetScalar(IN_MAXATOMS)+1e-2); - } - profile = (int)(mxGetScalar(IN_PROFILE)+1e-2); - - - /* check sizes */ - - if (D && x) { - n = mxGetM(IN_D); - m = mxGetN(IN_D); - L = mxGetN(IN_X); - - if (mxGetM(IN_X) != n) { - mexErrMsgTxt("D and X have incompatible sizes."); - } - - if (G) { - if (mxGetN(IN_G)!=mxGetM(IN_G)) { - mexErrMsgTxt("G must be a square matrix."); - } - if (mxGetN(IN_G) != m) { - mexErrMsgTxt("D and G have incompatible sizes."); - } - } - } - - else if (DtX && XtX) { - m = mxGetM(IN_DtX); - L = mxGetN(IN_DtX); - - /* set n to an arbitrary value that is at least the max possible number of selected atoms */ - - if (maxatoms>0) { - n = maxatoms; - } - else { - n = m; - } - - if ( !(mxGetM(IN_XtX)==L && mxGetN(IN_XtX)==1) && !(mxGetM(IN_XtX)==1 && mxGetN(IN_XtX)==L) ) { - mexErrMsgTxt("DtX and XtX have incompatible sizes."); - } - - if (mxGetN(IN_G)!=mxGetM(IN_G)) { - mexErrMsgTxt("G must be a square matrix."); - } - if (mxGetN(IN_G) != m) { - mexErrMsgTxt("DtX and G have incompatible sizes."); - } - } - - else { - mexErrMsgTxt("Either D and X, or DtX and XtX, must be specified."); - } - - - /* Do OMP! */ - - GAMMA_OUT = ompcoreGabor(D, x, DtX, XtX, G, n, m, L, maxatoms, eps, gmode, profile, msgdelta, 1); - - return; -}
--- a/util/ksvd utils/ompbox utils/omp2mexGabor.m Thu Jul 21 16:37:14 2011 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,23 +0,0 @@ -%This is the Matlab interface to the OMP2 MEX implementation. -%The function is not for independent use, only through omp2.m. - - -%OMP2MEX Matlab interface to the OMP2 MEX implementation. -% GAMMA = OMP2MEXGABOR(D,X,DtX,XtX,G,EPSILON,SPARSE_G,MSGDELTA,MAXATOMS,PROFILE) -% invokes the OMP2 MEX function according to the specified parameters. Not -% all the parameters are required. Those among D, X, DtX, XtX and G which -% are not specified should be passed as []. -% -% EPSILON - the target error. -% SPARSE_G - returns a sparse GAMMA when nonzero, full GAMMA when zero. -% MSGDELTA - the delay in secs between messages. Zero means no messages. -% MAXATOMS - the max number of atoms per signal, negative for no max. -% PROFILE - nonzero means that profiling information should be printed. - - -% Ron Rubinstein -% Computer Science Department -% Technion, Haifa 32000 Israel -% ronrubin@cs -% -% April 2009
--- a/util/ksvd utils/ompbox utils/ompGabor.m Thu Jul 21 16:37:14 2011 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,180 +0,0 @@ -function gamma = omp(varargin) -%OMP Sparsity-constrained Orthogonal Matching Pursuit. -% GAMMA = OMP(D,X,G,T) solves the optimization problem -% -% min |X - D*GAMMA|_2 s.t. |GAMMA|_0 <= T -% gamma -% -% for each of the signals in X, using Batch Orthogonal Matching Pursuit. -% Here, D is a dictionary with normalized columns, X is a matrix -% containing column signals, T is the # of non-zeros in each signal -% representation, and G is the Gramm matrix D'*D. The output GAMMA is a -% matrix containing the sparse representations as its columns. -% -% GAMMA = OMP(D,X,[],T) performs the same operation, but without the -% matrix G, using OMP-Cholesky. This call produces the same output as -% Batch-OMP, but is significantly slower. Using this syntax is only -% recommended when available memory is too small to store G. -% -% GAMMA = OMP(DtX,G,T) is the fastest implementation of OMP, but also -% requires the most memory. Here, DtX stores the projections D'*X. In this -% case Batch-OMP is used, but without having to compute D'*X in advance, -% which slightly improves runtime. Note that in general, the call -% -% GAMMA = OMP(D'*X,G,T); -% -% will be faster than the call -% -% GAMMA = OMP(D,X,G,T); -% -% due to optimized matrix multiplications in Matlab. However, when the -% entire matrix D'*X cannot be stored in memory, one of the other two -% versions can be used. Both compute D'*X for just one signal at a time, -% and thus require much less memory. -% -% GAMMA = OMP(...,PARAM1,VAL1,PARAM2,VAL2,...) specifies additional -% parameters for OMP. Available parameters are: -% -% 'gammamode' - Specifies the representation mode for GAMMA. Can be -% either 'full' or 'sparse', corresponding to a full or -% sparse matrix, respectively. By default, GAMMA is -% returned as a sparse matrix. -% 'messages' - Specifies whether progress messages should be displayed. -% When positive, this is the number of seconds between -% status prints. When negative, indicates that no messages -% should be displayed (this is the default). -% 'checkdict' - Specifies whether dictionary normalization should be -% verified. When set to 'on' (default) the dictionary -% atoms are verified to be of unit L2-norm. Setting this -% parameter to 'off' disables verification and accelerates -% function performance. Note that an unnormalized -% dictionary will produce invalid results. -% 'profile' - Can be either 'on' or 'off'. When 'on', profiling -% information is displayed at the end of the funciton -% execution. -% -% -% Summary of OMP versions: -% -% version | speed | memory -% -------------------------------------------------- -% OMP(DtX,G,T) | very fast | very large -% OMP(D,X,G,T) | fast | moderate -% OMP(D,X,[],T) | very slow | small -% -------------------------------------------------- -% -% -% References: -% [1] M. Elad, R. Rubinstein, and M. Zibulevsky, "Efficient Implementation -% of the K-SVD Algorithm using Batch Orthogonal Matching Pursuit", -% Technical Report - CS, Technion, April 2008. -% -% See also OMP2. - - -% Ron Rubinstein -% Computer Science Department -% Technion, Haifa 32000 Israel -% ronrubin@cs -% -% April 2009 - - -% default options - -sparse_gamma = 1; -msgdelta = -1; -checkdict = 1; -profile = 0; - - -% determine number of parameters - -paramnum = 1; -while (paramnum<=nargin && ~ischar(varargin{paramnum})) - paramnum = paramnum+1; -end -paramnum = paramnum-1; - - -% parse options - -for i = paramnum+1:2:length(varargin) - paramname = varargin{i}; - paramval = varargin{i+1}; - - switch lower(paramname) - - case 'gammamode' - if (strcmpi(paramval,'sparse')) - sparse_gamma = 1; - elseif (strcmpi(paramval,'full')) - sparse_gamma = 0; - else - error('Invalid GAMMA mode'); - end - - case 'messages' - msgdelta = paramval; - - case 'checkdict' - if (strcmpi(paramval,'on')) - checkdict = 1; - elseif (strcmpi(paramval,'off')) - checkdict = 0; - else - error('Invalid checkdict option'); - end - - case 'profile' - if (strcmpi(paramval,'on')) - profile = 1; - elseif (strcmpi(paramval,'off')) - profile = 0; - else - error('Invalid profile mode'); - end - - otherwise - error(['Unknown option: ' paramname]); - end - -end - - -% determine call type - -if (paramnum==3) - DtX = varargin{1}; - G = varargin{2}; - T = varargin{3}; - D = []; - X = []; -elseif (paramnum==4) - D = varargin{1}; - X = varargin{2}; - G = varargin{3}; - T = varargin{4}; - DtX = []; -else - error('Invalid number of parameters'); -end - - -% verify dictionary normalization - -if (checkdict) - if (isempty(G)) - atomnorms = sum(D.*D); - else - atomnorms = diag(G); - end - if (any(abs(atomnorms-1) > 1e-2)) - error('Dictionary columns must be normalized to unit length'); - end -end - - -% omp - -gamma = ompmexGabor(D,X,DtX,G,T,sparse_gamma,msgdelta,profile);
--- a/util/ksvd utils/ompbox utils/ompcore.c Thu Jul 21 16:37:14 2011 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,409 +0,0 @@ -/************************************************************************** - * - * File name: ompcore.c - * - * Ron Rubinstein - * Computer Science Department - * Technion, Haifa 32000 Israel - * ronrubin@cs - * - * Last Updated: 25.8.2009 - * - *************************************************************************/ - - -#include "ompcore.h" -#include "omputils.h" -#include "ompprof.h" -#include "myblas.h" -#include <math.h> -#include <string.h> - - - -/****************************************************************************** - * * - * Batch-OMP Implementation * - * * - ******************************************************************************/ - -mxArray* ompcore(double D[], double x[], double DtX[], double XtX[], double G[], mwSize n, mwSize m, mwSize L, - int T, double eps, int gamma_mode, int profile, double msg_delta, int erroromp) -{ - - profdata pd; - mxArray *Gamma; - mwIndex i, j, signum, pos, *ind, *gammaIr, *gammaJc, gamma_count; - mwSize allocated_coefs, allocated_cols; - int DtX_specified, XtX_specified, batchomp, standardomp, *selected_atoms; - double *alpha, *r, *Lchol, *c, *Gsub, *Dsub, sum, *gammaPr, *tempvec1, *tempvec2; - double eps2, resnorm, delta, deltaprev, secs_remain; - int mins_remain, hrs_remain; - clock_t lastprint_time, starttime; - - - - /*** status flags ***/ - - DtX_specified = (DtX!=0); /* indicates whether D'*x was provided */ - XtX_specified = (XtX!=0); /* indicates whether sum(x.*x) was provided */ - - standardomp = (G==0); /* batch-omp or standard omp are selected depending on availability of G */ - batchomp = !standardomp; - - - - /*** allocate output matrix ***/ - - - if (gamma_mode == FULL_GAMMA) { - - /* allocate full matrix of size m X L */ - - Gamma = mxCreateDoubleMatrix(m, L, mxREAL); - gammaPr = mxGetPr(Gamma); - gammaIr = 0; - gammaJc = 0; - } - else { - - /* allocate sparse matrix with room for allocated_coefs nonzeros */ - - /* for error-omp, begin with L*sqrt(n)/2 allocated nonzeros, otherwise allocate L*T nonzeros */ - allocated_coefs = erroromp ? (mwSize)(ceil(L*sqrt((double)n)/2.0) + 1.01) : L*T; - Gamma = mxCreateSparse(m, L, allocated_coefs, mxREAL); - gammaPr = mxGetPr(Gamma); - gammaIr = mxGetIr(Gamma); - gammaJc = mxGetJc(Gamma); - gamma_count = 0; - gammaJc[0] = 0; - } - - - /*** helper arrays ***/ - - alpha = (double*)mxMalloc(m*sizeof(double)); /* contains D'*residual */ - ind = (mwIndex*)mxMalloc(n*sizeof(mwIndex)); /* indices of selected atoms */ - selected_atoms = (int*)mxMalloc(m*sizeof(int)); /* binary array with 1's for selected atoms */ - c = (double*)mxMalloc(n*sizeof(double)); /* orthogonal projection result */ - - /* current number of columns in Dsub / Gsub / Lchol */ - allocated_cols = erroromp ? (mwSize)(ceil(sqrt((double)n)/2.0) + 1.01) : T; - - /* Cholesky decomposition of D_I'*D_I */ - Lchol = (double*)mxMalloc(n*allocated_cols*sizeof(double)); - - /* temporary vectors for various computations */ - tempvec1 = (double*)mxMalloc(m*sizeof(double)); - tempvec2 = (double*)mxMalloc(m*sizeof(double)); - - if (batchomp) { - /* matrix containing G(:,ind) - the columns of G corresponding to the selected atoms, in order of selection */ - Gsub = (double*)mxMalloc(m*allocated_cols*sizeof(double)); - } - else { - /* matrix containing D(:,ind) - the selected atoms from D, in order of selection */ - Dsub = (double*)mxMalloc(n*allocated_cols*sizeof(double)); - - /* stores the residual */ - r = (double*)mxMalloc(n*sizeof(double)); - } - - if (!DtX_specified) { - /* contains D'*x for the current signal */ - DtX = (double*)mxMalloc(m*sizeof(double)); - } - - - - /*** initializations for error omp ***/ - - if (erroromp) { - eps2 = eps*eps; /* compute eps^2 */ - if (T<0 || T>n) { /* unspecified max atom num - set max atoms to n */ - T = n; - } - } - - - - /*** initialize timers ***/ - - initprofdata(&pd); /* initialize profiling counters */ - starttime = clock(); /* record starting time for eta computations */ - lastprint_time = starttime; /* time of last status display */ - - - - /********************** perform omp for each signal **********************/ - - - - for (signum=0; signum<L; ++signum) { - - - /* initialize residual norm and deltaprev for error-omp */ - - if (erroromp) { - if (XtX_specified) { - resnorm = XtX[signum]; - } - else { - resnorm = dotprod(x+n*signum, x+n*signum, n); - addproftime(&pd, XtX_TIME); - } - deltaprev = 0; /* delta tracks the value of gamma'*G*gamma */ - } - else { - /* ignore residual norm stopping criterion */ - eps2 = 0; - resnorm = 1; - } - - - if (resnorm>eps2 && T>0) { - - /* compute DtX */ - - if (!DtX_specified) { - matT_vec(1, D, x+n*signum, DtX, n, m); - addproftime(&pd, DtX_TIME); - } - - - /* initialize alpha := DtX */ - - memcpy(alpha, DtX + m*signum*DtX_specified, m*sizeof(double)); - - - /* mark all atoms as unselected */ - - for (i=0; i<m; ++i) { - selected_atoms[i] = 0; - } - - } - - - /* main loop */ - - i=0; - while (resnorm>eps2 && i<T) { - - /* index of next atom */ - - pos = maxabs(alpha, m); - addproftime(&pd, MAXABS_TIME); - - - /* stop criterion: selected same atom twice, or inner product too small */ - - if (selected_atoms[pos] || alpha[pos]*alpha[pos]<1e-14) { - break; - } - - - /* mark selected atom */ - - ind[i] = pos; - selected_atoms[pos] = 1; - - - /* matrix reallocation */ - - if (erroromp && i>=allocated_cols) { - - allocated_cols = (mwSize)(ceil(allocated_cols*MAT_INC_FACTOR) + 1.01); - - Lchol = (double*)mxRealloc(Lchol,n*allocated_cols*sizeof(double)); - - batchomp ? (Gsub = (double*)mxRealloc(Gsub,m*allocated_cols*sizeof(double))) : - (Dsub = (double*)mxRealloc(Dsub,n*allocated_cols*sizeof(double))) ; - } - - - /* append column to Gsub or Dsub */ - - if (batchomp) { - memcpy(Gsub+i*m, G+pos*m, m*sizeof(double)); - } - else { - memcpy(Dsub+i*n, D+pos*n, n*sizeof(double)); - } - - - /*** Cholesky update ***/ - - if (i==0) { - *Lchol = 1; - } - else { - - /* incremental Cholesky decomposition: compute next row of Lchol */ - - if (standardomp) { - matT_vec(1, Dsub, D+n*pos, tempvec1, n, i); /* compute tempvec1 := Dsub'*d where d is new atom */ - addproftime(&pd, DtD_TIME); - } - else { - vec_assign(tempvec1, Gsub+i*m, ind, i); /* extract tempvec1 := Gsub(ind,i) */ - } - backsubst('L', Lchol, tempvec1, tempvec2, n, i); /* compute tempvec2 = Lchol \ tempvec1 */ - for (j=0; j<i; ++j) { /* write tempvec2 to end of Lchol */ - Lchol[j*n+i] = tempvec2[j]; - } - - /* compute Lchol(i,i) */ - sum = 0; - for (j=0; j<i; ++j) { /* compute sum of squares of last row without Lchol(i,i) */ - sum += SQR(Lchol[j*n+i]); - } - if ( (1-sum) <= 1e-14 ) { /* Lchol(i,i) is zero => selected atoms are dependent */ - break; - } - Lchol[i*n+i] = sqrt(1-sum); - } - - addproftime(&pd, LCHOL_TIME); - - i++; - - - /* perform orthogonal projection and compute sparse coefficients */ - - vec_assign(tempvec1, DtX + m*signum*DtX_specified, ind, i); /* extract tempvec1 = DtX(ind) */ - cholsolve('L', Lchol, tempvec1, c, n, i); /* solve LL'c = tempvec1 for c */ - addproftime(&pd, COMPCOEF_TIME); - - - /* update alpha = D'*residual */ - - if (standardomp) { - mat_vec(-1, Dsub, c, r, n, i); /* compute r := -Dsub*c */ - vec_sum(1, x+n*signum, r, n); /* compute r := x+r */ - - - /*memcpy(r, x+n*signum, n*sizeof(double)); /* assign r := x */ - /*mat_vec1(-1, Dsub, c, 1, r, n, i); /* compute r := r-Dsub*c */ - - addproftime(&pd, COMPRES_TIME); - matT_vec(1, D, r, alpha, n, m); /* compute alpha := D'*r */ - addproftime(&pd, DtR_TIME); - - /* update residual norm */ - if (erroromp) { - resnorm = dotprod(r, r, n); - addproftime(&pd, UPDATE_RESNORM_TIME); - } - } - else { - mat_vec(1, Gsub, c, tempvec1, m, i); /* compute tempvec1 := Gsub*c */ - memcpy(alpha, DtX + m*signum*DtX_specified, m*sizeof(double)); /* set alpha = D'*x */ - vec_sum(-1, tempvec1, alpha, m); /* compute alpha := alpha - tempvec1 */ - addproftime(&pd, UPDATE_DtR_TIME); - - /* update residual norm */ - if (erroromp) { - vec_assign(tempvec2, tempvec1, ind, i); /* assign tempvec2 := tempvec1(ind) */ - delta = dotprod(c,tempvec2,i); /* compute c'*tempvec2 */ - resnorm = resnorm - delta + deltaprev; /* residual norm update */ - deltaprev = delta; - addproftime(&pd, UPDATE_RESNORM_TIME); - } - } - } - - - /*** generate output vector gamma ***/ - - if (gamma_mode == FULL_GAMMA) { /* write the coefs in c to their correct positions in gamma */ - for (j=0; j<i; ++j) { - gammaPr[m*signum + ind[j]] = c[j]; - } - } - else { - /* sort the coefs by index before writing them to gamma */ - quicksort(ind,c,i); - addproftime(&pd, INDEXSORT_TIME); - - /* gamma is full - reallocate */ - if (gamma_count+i >= allocated_coefs) { - - while(gamma_count+i >= allocated_coefs) { - allocated_coefs = (mwSize)(ceil(GAMMA_INC_FACTOR*allocated_coefs) + 1.01); - } - - mxSetNzmax(Gamma, allocated_coefs); - mxSetPr(Gamma, mxRealloc(gammaPr, allocated_coefs*sizeof(double))); - mxSetIr(Gamma, mxRealloc(gammaIr, allocated_coefs*sizeof(mwIndex))); - - gammaPr = mxGetPr(Gamma); - gammaIr = mxGetIr(Gamma); - } - - /* append coefs to gamma and update the indices */ - for (j=0; j<i; ++j) { - gammaPr[gamma_count] = c[j]; - gammaIr[gamma_count] = ind[j]; - gamma_count++; - } - gammaJc[signum+1] = gammaJc[signum] + i; - } - - - - /*** display status messages ***/ - - if (msg_delta>0 && (clock()-lastprint_time)/(double)CLOCKS_PER_SEC >= msg_delta) - { - lastprint_time = clock(); - - /* estimated remainig time */ - secs2hms( ((L-signum-1)/(double)(signum+1)) * ((lastprint_time-starttime)/(double)CLOCKS_PER_SEC) , - &hrs_remain, &mins_remain, &secs_remain); - - mexPrintf("omp: signal %d / %d, estimated remaining time: %02d:%02d:%05.2f\n", - signum+1, L, hrs_remain, mins_remain, secs_remain); - mexEvalString("drawnow;"); - } - - } - - /* end omp */ - - - - /*** print final messages ***/ - - if (msg_delta>0) { - mexPrintf("omp: signal %d / %d\n", signum, L); - } - - if (profile) { - printprofinfo(&pd, erroromp, batchomp, L); - } - - - - /* free memory */ - - if (!DtX_specified) { - mxFree(DtX); - } - if (standardomp) { - mxFree(r); - mxFree(Dsub); - } - else { - mxFree(Gsub); - } - mxFree(tempvec2); - mxFree(tempvec1); - mxFree(Lchol); - mxFree(c); - mxFree(selected_atoms); - mxFree(ind); - mxFree(alpha); - - return Gamma; -}
--- a/util/ksvd utils/ompbox utils/ompcore.h Thu Jul 21 16:37:14 2011 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,80 +0,0 @@ -/************************************************************************** - * - * File name: ompcore.h - * - * Ron Rubinstein - * Computer Science Department - * Technion, Haifa 32000 Israel - * ronrubin@cs - * - * Last Updated: 18.8.2009 - * - * Contains the core implementation of Batch-OMP / OMP-Cholesky. - * - *************************************************************************/ - - -#ifndef __OMP_CORE_H__ -#define __OMP_CORE_H__ - - -#include "mex.h" - - - -/************************************************************************** - * Perform Batch-OMP or OMP-Cholesky on a specified set of signals, using - * either a fixed number of atoms or an error bound. - * - * Parameters (not all required): - * - * D - the dictionary, of size n X m - * x - the signals, of size n X L - * DtX - D'*x, of size m X L - * XtX - squared norms of the signals in x, sum(x.*x), of length L - * G - D'*D, of size m X m - * T - target sparsity, or maximal number of atoms for error-based OMP - * eps - target residual norm for error-based OMP - * gamma_mode - one of the constants FULL_GAMMA or SPARSE_GAMMA - * profile - if non-zero, profiling info is printed - * msg_delta - positive: the # of seconds between status prints, otherwise: nothing is printed - * erroromp - if nonzero indicates error-based OMP, otherwise fixed sparsity OMP - * - * Usage: - * - * The function can be called using different parameters, and will have - * different complexity depending on the parameters specified. Arrays which - * are not specified should be passed as null (0). When G is specified, - * Batch-OMP is performed. Otherwise, OMP-Cholesky is performed. - * - * Fixed-sparsity usage: - * --------------------- - * Either DtX, or D and x, must be specified. Specifying DtX is more efficient. - * XtX does not need to be specified. - * When D and x are specified, G is not required. However, not providing G - * will significantly degrade efficiency. - * The number of atoms must be specified in T. The value of eps is ignored. - * Finally, set erroromp to 0. - * - * Error-OMP usage: - * ---------------- - * Either DtX and Xtx, or D and x, must be specified. Specifying DtX and XtX - * is more efficient. - * When D and x are specified, G is not required. However, not providing G - * will significantly degrade efficiency. - * The target error must be specified in eps. A hard limit on the number - * of atoms can also be specified via the parameter T. Otherwise, T should - * be negative. Finally, set erroromp to nonzero. - * - * - * Returns: - * An mxArray containing the sparse representations of the signals in x - * (allocated using the appropriate mxCreateXXX() function). - * The array is either full or sparse, depending on gamma_mode. - * - **************************************************************************/ -mxArray* ompcore(double D[], double x[], double DtX[], double XtX[], double G[], mwSize n, mwSize m, mwSize L, - int T, double eps, int gamma_mode, int profile, double msg_delta, int erroromp); - - -#endif
--- a/util/ksvd utils/ompbox utils/ompcoreGabor.c Thu Jul 21 16:37:14 2011 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,465 +0,0 @@ -/************************************************************************** - * - * File name: ompcoreGabor.c - * - * Ron Rubinstein - * Computer Science Department - * Technion, Haifa 32000 Israel - * ronrubin@cs - * - * Last Updated: 25.8.2009 - * - * Modified by Ivan damnjanovic July 2011 - * Takes to atoms per iteration. It should be used for Gabor dictionaries - * as specified in - * "Audio Inpainting" Amir Adler, Valentin Emiya, Maria G. Jafari, - * Michael Elad, Remi Gribonval and Mark D. Plumbley - * Draft version: March 6, 2011 - * - *************************************************************************/ - - -#include "ompcoreGabor.h" -#include "omputils.h" -#include "ompprof.h" -#include "myblas.h" -#include <math.h> -#include <string.h> - - - -/****************************************************************************** - * * - * Batch-OMP Implementation * - * * - ******************************************************************************/ - -mxArray* ompcoreGabor(double D[], double x[], double DtX[], double XtX[], double G[], mwSize n, mwSize m, mwSize L, - int T, double eps, int gamma_mode, int profile, double msg_delta, int erroromp) -{ - - profdata pd; - mxArray *Gamma; - mwIndex i, j, k, signum, pos, *ind, *gammaIr, *gammaJc, gamma_count; - mwSize allocated_coefs, allocated_cols; - int DtX_specified, XtX_specified, batchomp, standardomp, *selected_atoms; - double *proj, *proj1, *proj2, *D1, *D2, *D1D2, *n12, *alpha, *beta, *error; - double *r, *Lchol, *c, *Gsub, *Dsub, sum, *gammaPr, *tempvec1, *tempvec2; - double eps2, resnorm, delta, deltaprev, secs_remain; - int mins_remain, hrs_remain; - clock_t lastprint_time, starttime; - - - - /*** status flags ***/ - - DtX_specified = (DtX!=0); /* indicates whether D'*x was provided */ - XtX_specified = (XtX!=0); /* indicates whether sum(x.*x) was provided */ - - standardomp = 1;//(G==0); /* batch-omp or standard omp are selected depending on availability of G */ - batchomp = !standardomp; - - - - /*** allocate output matrix ***/ - - - if (gamma_mode == FULL_GAMMA) { - - /* allocate full matrix of size m X L */ - - Gamma = mxCreateDoubleMatrix(m, L, mxREAL); - gammaPr = mxGetPr(Gamma); - gammaIr = 0; - gammaJc = 0; - } - else { - - /* allocate sparse matrix with room for allocated_coefs nonzeros */ - - /* for error-omp, begin with L*sqrt(n)/2 allocated nonzeros, otherwise allocate L*T nonzeros */ - allocated_coefs = erroromp ? (mwSize)(ceil(L*sqrt((double)n)/2.0) + 1.01) : L*T; - Gamma = mxCreateSparse(m, L, allocated_coefs, mxREAL); - gammaPr = mxGetPr(Gamma); - gammaIr = mxGetIr(Gamma); - gammaJc = mxGetJc(Gamma); - gamma_count = 0; - gammaJc[0] = 0; - } - - - /*** helper arrays ***/ - /* Ivan Damnjanovic July 2011*/ - proj = (double*)mxMalloc(m*sizeof(double)); - proj1 = (double*)mxMalloc(m/2*sizeof(double)); - proj2 = (double*)mxMalloc(m/2*sizeof(double)); - D1 = (double*)mxMalloc(n*m/2*sizeof(double)); - D2 = (double*)mxMalloc(n*m/2*sizeof(double)); - memcpy(D1, D , n*m/2*sizeof(double)); - memcpy(D2, D+n*m/2, n*m/2*sizeof(double)); - - D1D2 = (double*)mxMalloc(m/2*sizeof(double)); - n12 = (double*)mxMalloc(m/2*sizeof(double)); - - vec_smult(1,D2, D1, n*m/2); - - for (i=0; i<m/2; i++) { - D1D2[i]=0; - n12[i]=0; - for (j=0; j<n; j++) { - D1D2[i] += D1[i*n+j]; - } - n12[i]=1/(1-D1D2[i]*D1D2[i]); - } - - memcpy(D1, D , n*m/2*sizeof(double)); - - alpha = (double*)mxMalloc(m/2*sizeof(double)); /* contains D'*residual */ - beta = (double*)mxMalloc(m/2*sizeof(double)); - error = (double*)mxMalloc(m/2*sizeof(double)); - - ind = (mwIndex*)mxMalloc(m*sizeof(mwIndex)); /* indices of selected atoms */ - selected_atoms = (int*)mxMalloc(m*sizeof(int)); /* binary array with 1's for selected atoms */ - c = (double*)mxMalloc(n*sizeof(double)); /* orthogonal projection result */ - - /* current number of columns in Dsub / Gsub / Lchol */ - allocated_cols = erroromp ? (mwSize)(ceil(sqrt((double)n)/2.0) + 1.01) : T; - - /* Cholesky decomposition of D_I'*D_I */ - Lchol = (double*)mxMalloc(n*allocated_cols*sizeof(double)); - - /* temporary vectors for various computations */ - tempvec1 = (double*)mxMalloc(m*sizeof(double)); - tempvec2 = (double*)mxMalloc(m*sizeof(double)); - - if (batchomp) { - /* matrix containing G(:,ind) - the columns of G corresponding to the selected atoms, in order of selection */ - Gsub = (double*)mxMalloc(m*allocated_cols*sizeof(double)); - } - else { - /* matrix containing D(:,ind) - the selected atoms from D, in order of selection */ - Dsub = (double*)mxMalloc(n*allocated_cols*sizeof(double)); - - /* stores the residual */ - r = (double*)mxMalloc(n*sizeof(double)); - } - - if (!DtX_specified) { - /* contains D'*x for the current signal */ - DtX = (double*)mxMalloc(m*sizeof(double)); - } - - - - /*** initializations for error omp ***/ - - if (erroromp) { - eps2 = eps*eps; /* compute eps^2 */ - if (T<0 || T>n) { /* unspecified max atom num - set max atoms to n */ - T = n; - } - } - - - - /*** initialize timers ***/ - - initprofdata(&pd); /* initialize profiling counters */ - starttime = clock(); /* record starting time for eta computations */ - lastprint_time = starttime; /* time of last status display */ - - - - /********************** perform omp for each signal **********************/ - - - - for (signum=0; signum<L; ++signum) { - - - /* initialize residual norm and deltaprev for error-omp */ - - if (erroromp) { - if (XtX_specified) { - resnorm = XtX[signum]; - } - else { - resnorm = dotprod(x+n*signum, x+n*signum, n); - addproftime(&pd, XtX_TIME); - } - deltaprev = 0; /* delta tracks the value of gamma'*G*gamma */ - } - else { - /* ignore residual norm stopping criterion */ - eps2 = 0; - resnorm = 1; - } - - - if (resnorm>eps2 && T>0) { - - /* compute DtX */ - - if (!DtX_specified) { - matT_vec(1, D, x+n*signum, DtX, n, m); - addproftime(&pd, DtX_TIME); - memcpy(r , x+n*signum, n*sizeof(double)); - } - - - /* initialize projections to D1 and D2 := DtX */ - - memcpy(proj, DtX + m*signum*DtX_specified, m*sizeof(double)); - - - /* mark all atoms as unselected */ - - for (i=0; i<m; ++i) { - selected_atoms[i] = 0; - } - - } - - - /* main loop */ - - i=0; - while (resnorm>eps2 && i<T) { - - /* index of next atom */ - memcpy(proj1, proj, m/2*sizeof(double)); - memcpy(proj2, proj + m/2, m/2*sizeof(double)); - for (k=0; k<m/2; k++){ - alpha[k] = (proj1[k] - D1D2[k]*proj2[k])*n12[k]; - beta[k] = (proj2[k] - D1D2[k]*proj1[k])*n12[k]; - } - for (k=0; k<m/2; k++){ - error[k]=0; - for (j=0; j<n; j++){ - error[k]+= (abs(r[j])-D1[k*n+j]*alpha[k]-D2[k*n+j]*beta[k])*(abs(r[j])-D1[k*n+j]*alpha[k]-D2[k*n+j]*beta[k]); - } - } - pos = maxabs(error, m/2); - addproftime(&pd, MAXABS_TIME); - - - /* stop criterion: selected same atom twice, or inner product too small */ - - if (selected_atoms[pos] || alpha[pos]*alpha[pos]<1e-14) { - break; - } - - for (k=0;k<2;k++){ - /* mark selected atom */ - - ind[i] = pos+k*m/2; - selected_atoms[pos+k*m/2] = 1; - - - /* matrix reallocation */ - - if (erroromp && i>=allocated_cols) { - - allocated_cols = (mwSize)(ceil(allocated_cols*MAT_INC_FACTOR) + 1.01); - - Lchol = (double*)mxRealloc(Lchol,n*allocated_cols*sizeof(double)); - - batchomp ? (Gsub = (double*)mxRealloc(Gsub,m*allocated_cols*sizeof(double))) : - (Dsub = (double*)mxRealloc(Dsub,n*allocated_cols*sizeof(double))) ; - } - - - /* append column to Gsub or Dsub */ - - if (batchomp) { - memcpy(Gsub+i*m, G+(pos+k*m/2)*m, m*sizeof(double)); - } - else { - memcpy(Dsub+(i)*n, D+(pos+k*m/2)*n, n*sizeof(double)); - } - - - /*** Cholesky update ***/ - - if (i==0) { - *Lchol = 1; - } - else { - - /* incremental Cholesky decomposition: compute next row of Lchol */ - - if (standardomp) { - matT_vec(1, Dsub, D+n*(pos+k*m/2), tempvec1, n, i); /* compute tempvec1 := Dsub'*d where d is new atom */ - addproftime(&pd, DtD_TIME); - } - else { - vec_assign(tempvec1, Gsub+i*m, ind, i); /* extract tempvec1 := Gsub(ind,i) */ - } - backsubst('L', Lchol, tempvec1, tempvec2, n, i); /* compute tempvec2 = Lchol \ tempvec1 */ - for (j=0; j<i; ++j) { /* write tempvec2 to end of Lchol */ - Lchol[j*n+i] = tempvec2[j]; - } - - /* compute Lchol(i,i) */ - sum = 0; - for (j=0; j<i; ++j) { /* compute sum of squares of last row without Lchol(i,i) */ - sum += SQR(Lchol[j*n+i]); - } - if ( (1-sum) <= 1e-14 ) { /* Lchol(i,i) is zero => selected atoms are dependent */ - break; - } - Lchol[i*n+i] = sqrt(1-sum); - } - - addproftime(&pd, LCHOL_TIME); - - i++; - - } - /* perform orthogonal projection and compute sparse coefficients */ - - vec_assign(tempvec1, DtX + m*signum*DtX_specified, ind, i); /* extract tempvec1 = DtX(ind) */ - cholsolve('L', Lchol, tempvec1, c, n, i); /* solve LL'c = tempvec1 for c */ - addproftime(&pd, COMPCOEF_TIME); - - - /* update alpha = D'*residual */ - - if (standardomp) { - mat_vec(-1, Dsub, c, r, n, i); /* compute r := -Dsub*c */ - vec_sum(1, x+n*signum, r, n); /* compute r := x+r */ - - - /*memcpy(r, x+n*signum, n*sizeof(double)); /* assign r := x */ - /*mat_vec1(-1, Dsub, c, 1, r, n, i); /* compute r := r-Dsub*c */ - - addproftime(&pd, COMPRES_TIME); - matT_vec(1, D, r, proj, n, m); /* compute proj := D'*r */ - addproftime(&pd, DtR_TIME); - - /* update residual norm */ - if (erroromp) { - resnorm = dotprod(r, r, n); - addproftime(&pd, UPDATE_RESNORM_TIME); - } - } - else { - mat_vec(1, Gsub, c, tempvec1, m, i); /* compute tempvec1 := Gsub*c */ - memcpy(proj, DtX + m*signum*DtX_specified, m*sizeof(double)); /* set proj = D'*x */ - vec_sum(-1, tempvec1, proj, m); /* compute proj := proj - tempvec1 */ - addproftime(&pd, UPDATE_DtR_TIME); - - /* update residual norm */ - if (erroromp) { - vec_assign(tempvec2, tempvec1, ind, i); /* assign tempvec2 := tempvec1(ind) */ - delta = dotprod(c,tempvec2,i); /* compute c'*tempvec2 */ - resnorm = resnorm - delta + deltaprev; /* residual norm update */ - deltaprev = delta; - addproftime(&pd, UPDATE_RESNORM_TIME); - } - } - } - - - /*** generate output vector gamma ***/ - - if (gamma_mode == FULL_GAMMA) { /* write the coefs in c to their correct positions in gamma */ - for (j=0; j<i; ++j) { - gammaPr[m*signum + ind[j]] = c[j]; - } - } - else { - /* sort the coefs by index before writing them to gamma */ - quicksort(ind,c,i); - addproftime(&pd, INDEXSORT_TIME); - - /* gamma is full - reallocate */ - if (gamma_count+i >= allocated_coefs) { - - while(gamma_count+i >= allocated_coefs) { - allocated_coefs = (mwSize)(ceil(GAMMA_INC_FACTOR*allocated_coefs) + 1.01); - } - - mxSetNzmax(Gamma, allocated_coefs); - mxSetPr(Gamma, mxRealloc(gammaPr, allocated_coefs*sizeof(double))); - mxSetIr(Gamma, mxRealloc(gammaIr, allocated_coefs*sizeof(mwIndex))); - - gammaPr = mxGetPr(Gamma); - gammaIr = mxGetIr(Gamma); - } - - /* append coefs to gamma and update the indices */ - for (j=0; j<i; ++j) { - gammaPr[gamma_count] = c[j]; - gammaIr[gamma_count] = ind[j]; - gamma_count++; - } - gammaJc[signum+1] = gammaJc[signum] + i; - } - - - - /*** display status messages ***/ - - if (msg_delta>0 && (clock()-lastprint_time)/(double)CLOCKS_PER_SEC >= msg_delta) - { - lastprint_time = clock(); - - /* estimated remainig time */ - secs2hms( ((L-signum-1)/(double)(signum+1)) * ((lastprint_time-starttime)/(double)CLOCKS_PER_SEC) , - &hrs_remain, &mins_remain, &secs_remain); - - mexPrintf("omp: signal %d / %d, estimated remaining time: %02d:%02d:%05.2f\n", - signum+1, L, hrs_remain, mins_remain, secs_remain); - mexEvalString("drawnow;"); - } - - } - - /* end omp */ - - - - /*** print final messages ***/ - - if (msg_delta>0) { - mexPrintf("omp: signal %d / %d\n", signum, L); - } - - if (profile) { - printprofinfo(&pd, erroromp, batchomp, L); - } - - - - /* free memory */ - - if (!DtX_specified) { - mxFree(DtX); - } - if (standardomp) { - mxFree(r); - mxFree(Dsub); - } - else { - mxFree(Gsub); - } - mxFree(tempvec2); - mxFree(tempvec1); - mxFree(Lchol); - mxFree(c); - mxFree(selected_atoms); - mxFree(ind); - mxFree(proj); - mxFree(proj1); - mxFree(proj2); - mxFree(D1); - mxFree(D2); - mxFree(D1D2); - mxFree(n12); - mxFree(alpha); - mxFree(beta); - mxFree(error); - - return Gamma; -}
--- a/util/ksvd utils/ompbox utils/ompcoreGabor.h Thu Jul 21 16:37:14 2011 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,80 +0,0 @@ -/************************************************************************** - * - * File name: ompcore.h - * - * Ron Rubinstein - * Computer Science Department - * Technion, Haifa 32000 Israel - * ronrubin@cs - * - * Last Updated: 18.8.2009 - * - * Contains the core implementation of Batch-OMP / OMP-Cholesky. - * - *************************************************************************/ - - -#ifndef __OMP_CORE_H__ -#define __OMP_CORE_H__ - - -#include "mex.h" - - - -/************************************************************************** - * Perform Batch-OMP or OMP-Cholesky on a specified set of signals, using - * either a fixed number of atoms or an error bound. - * - * Parameters (not all required): - * - * D - the dictionary, of size n X m - * x - the signals, of size n X L - * DtX - D'*x, of size m X L - * XtX - squared norms of the signals in x, sum(x.*x), of length L - * G - D'*D, of size m X m - * T - target sparsity, or maximal number of atoms for error-based OMP - * eps - target residual norm for error-based OMP - * gamma_mode - one of the constants FULL_GAMMA or SPARSE_GAMMA - * profile - if non-zero, profiling info is printed - * msg_delta - positive: the # of seconds between status prints, otherwise: nothing is printed - * erroromp - if nonzero indicates error-based OMP, otherwise fixed sparsity OMP - * - * Usage: - * - * The function can be called using different parameters, and will have - * different complexity depending on the parameters specified. Arrays which - * are not specified should be passed as null (0). When G is specified, - * Batch-OMP is performed. Otherwise, OMP-Cholesky is performed. - * - * Fixed-sparsity usage: - * --------------------- - * Either DtX, or D and x, must be specified. Specifying DtX is more efficient. - * XtX does not need to be specified. - * When D and x are specified, G is not required. However, not providing G - * will significantly degrade efficiency. - * The number of atoms must be specified in T. The value of eps is ignored. - * Finally, set erroromp to 0. - * - * Error-OMP usage: - * ---------------- - * Either DtX and Xtx, or D and x, must be specified. Specifying DtX and XtX - * is more efficient. - * When D and x are specified, G is not required. However, not providing G - * will significantly degrade efficiency. - * The target error must be specified in eps. A hard limit on the number - * of atoms can also be specified via the parameter T. Otherwise, T should - * be negative. Finally, set erroromp to nonzero. - * - * - * Returns: - * An mxArray containing the sparse representations of the signals in x - * (allocated using the appropriate mxCreateXXX() function). - * The array is either full or sparse, depending on gamma_mode. - * - **************************************************************************/ -mxArray* ompcoreGabor(double D[], double x[], double DtX[], double XtX[], double G[], mwSize n, mwSize m, mwSize L, - int T, double eps, int gamma_mode, int profile, double msg_delta, int erroromp); - - -#endif
--- a/util/ksvd utils/ompbox utils/ompmex.c Thu Jul 21 16:37:14 2011 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,133 +0,0 @@ -/************************************************************************** - * - * File name: ompmex.c - * - * Ron Rubinstein - * Computer Science Department - * Technion, Haifa 32000 Israel - * ronrubin@cs - * - * Last Updated: 18.8.2009 - * - *************************************************************************/ - -#include "ompcore.h" -#include "omputils.h" -#include "mexutils.h" - - -/* Input Arguments */ - -#define IN_D prhs[0] -#define IN_X prhs[1] -#define IN_DtX prhs[2] -#define IN_G prhs[3] -#define IN_T prhs[4] -#define IN_SPARSE_G prhs[5] -#define IN_MSGDELTA prhs[6] -#define IN_PROFILE prhs[7] - - -/* Output Arguments */ - -#define GAMMA_OUT plhs[0] - - -/***************************************************************************************/ - - -void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray*prhs[]) - -{ - double *D, *x, *DtX, *G, msgdelta; - int gmode, profile, T; - mwSize m, n, L; /* D is n x m , X is n x L, DtX is m x L */ - - - /* check parameters */ - - checkmatrix(IN_D, "OMP", "D"); - checkmatrix(IN_X, "OMP", "X"); - checkmatrix(IN_DtX, "OMP", "DtX"); - checkmatrix(IN_G, "OMP", "G"); - - checkscalar(IN_T, "OMP", "T"); - checkscalar(IN_SPARSE_G, "OMP", "sparse_g"); - checkscalar(IN_MSGDELTA, "OMP", "msgdelta"); - checkscalar(IN_PROFILE, "OMP", "profile"); - - - /* get parameters */ - - x = D = DtX = G = 0; - - if (!mxIsEmpty(IN_D)) - D = mxGetPr(IN_D); - - if (!mxIsEmpty(IN_X)) - x = mxGetPr(IN_X); - - if (!mxIsEmpty(IN_DtX)) - DtX = mxGetPr(IN_DtX); - - if (!mxIsEmpty(IN_G)) - G = mxGetPr(IN_G); - - T = (int)(mxGetScalar(IN_T)+1e-2); - if ((int)(mxGetScalar(IN_SPARSE_G)+1e-2)) { - gmode = SPARSE_GAMMA; - } - else { - gmode = FULL_GAMMA; - } - msgdelta = mxGetScalar(IN_MSGDELTA); - profile = (int)(mxGetScalar(IN_PROFILE)+1e-2); - - - /* check sizes */ - - if (D && x) { - n = mxGetM(IN_D); - m = mxGetN(IN_D); - L = mxGetN(IN_X); - - if (mxGetM(IN_X) != n) { - mexErrMsgTxt("D and X have incompatible sizes."); - } - - if (G) { - if (mxGetN(IN_G)!=mxGetM(IN_G)) { - mexErrMsgTxt("G must be a square matrix."); - } - if (mxGetN(IN_G) != m) { - mexErrMsgTxt("D and G have incompatible sizes."); - } - } - } - - else if (DtX) { - m = mxGetM(IN_DtX); - L = mxGetN(IN_DtX); - - n = T; /* arbitrary - it is enough to assume signal length is T */ - - if (mxGetN(IN_G)!=mxGetM(IN_G)) { - mexErrMsgTxt("G must be a square matrix."); - } - if (mxGetN(IN_G) != m) { - mexErrMsgTxt("DtX and G have incompatible sizes."); - } - } - - else { - mexErrMsgTxt("Either D and X, or DtX, must be specified."); - } - - - /* Do OMP! */ - - GAMMA_OUT = ompcore(D, x, DtX, 0, G, n, m, L, T, 0, gmode, profile, msgdelta, 0); - - return; -} -
--- a/util/ksvd utils/ompbox utils/ompmex.m Thu Jul 21 16:37:14 2011 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,22 +0,0 @@ -%This is the Matlab interface to the OMP MEX implementation. -%The function is not for independent use, only through omp.m. - - -%OMPMEX Matlab interface to the OMP MEX implementation. -% GAMMA = OMPMEX(D,X,DtX,G,L,SPARSE_G,MSGDELTA,PROFILE) invokes the OMP -% MEX function according to the specified parameters. Not all the -% parameters are required. Those among D, X, DtX and G which are not -% specified should be passed as []. -% -% L - the target sparsity. -% SPARSE_G - returns a sparse GAMMA when nonzero, full GAMMA when zero. -% MSGDELTA - the delay in secs between messages. Zero means no messages. -% PROFILE - nonzero means that profiling information should be printed. - - -% Ron Rubinstein -% Computer Science Department -% Technion, Haifa 32000 Israel -% ronrubin@cs -% -% April 2009
--- a/util/ksvd utils/ompbox utils/ompmexGabor.c Thu Jul 21 16:37:14 2011 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,133 +0,0 @@ -/************************************************************************** - * - * File name: ompmex.c - * - * Ron Rubinstein - * Computer Science Department - * Technion, Haifa 32000 Israel - * ronrubin@cs - * - * Last Updated: 18.8.2009 - * - *************************************************************************/ - -#include "ompcoreGabor.h" -#include "omputils.h" -#include "mexutils.h" - - -/* Input Arguments */ - -#define IN_D prhs[0] -#define IN_X prhs[1] -#define IN_DtX prhs[2] -#define IN_G prhs[3] -#define IN_T prhs[4] -#define IN_SPARSE_G prhs[5] -#define IN_MSGDELTA prhs[6] -#define IN_PROFILE prhs[7] - - -/* Output Arguments */ - -#define GAMMA_OUT plhs[0] - - -/***************************************************************************************/ - - -void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray*prhs[]) - -{ - double *D, *x, *DtX, *G, msgdelta; - int gmode, profile, T; - mwSize m, n, L; /* D is n x m , X is n x L, DtX is m x L */ - - - /* check parameters */ - - checkmatrix(IN_D, "OMP", "D"); - checkmatrix(IN_X, "OMP", "X"); - checkmatrix(IN_DtX, "OMP", "DtX"); - checkmatrix(IN_G, "OMP", "G"); - - checkscalar(IN_T, "OMP", "T"); - checkscalar(IN_SPARSE_G, "OMP", "sparse_g"); - checkscalar(IN_MSGDELTA, "OMP", "msgdelta"); - checkscalar(IN_PROFILE, "OMP", "profile"); - - - /* get parameters */ - - x = D = DtX = G = 0; - - if (!mxIsEmpty(IN_D)) - D = mxGetPr(IN_D); - - if (!mxIsEmpty(IN_X)) - x = mxGetPr(IN_X); - - if (!mxIsEmpty(IN_DtX)) - DtX = mxGetPr(IN_DtX); - - if (!mxIsEmpty(IN_G)) - G = mxGetPr(IN_G); - - T = (int)(mxGetScalar(IN_T)+1e-2); - if ((int)(mxGetScalar(IN_SPARSE_G)+1e-2)) { - gmode = SPARSE_GAMMA; - } - else { - gmode = FULL_GAMMA; - } - msgdelta = mxGetScalar(IN_MSGDELTA); - profile = (int)(mxGetScalar(IN_PROFILE)+1e-2); - - - /* check sizes */ - - if (D && x) { - n = mxGetM(IN_D); - m = mxGetN(IN_D); - L = mxGetN(IN_X); - - if (mxGetM(IN_X) != n) { - mexErrMsgTxt("D and X have incompatible sizes."); - } - - if (G) { - if (mxGetN(IN_G)!=mxGetM(IN_G)) { - mexErrMsgTxt("G must be a square matrix."); - } - if (mxGetN(IN_G) != m) { - mexErrMsgTxt("D and G have incompatible sizes."); - } - } - } - - else if (DtX) { - m = mxGetM(IN_DtX); - L = mxGetN(IN_DtX); - - n = T; /* arbitrary - it is enough to assume signal length is T */ - - if (mxGetN(IN_G)!=mxGetM(IN_G)) { - mexErrMsgTxt("G must be a square matrix."); - } - if (mxGetN(IN_G) != m) { - mexErrMsgTxt("DtX and G have incompatible sizes."); - } - } - - else { - mexErrMsgTxt("Either D and X, or DtX, must be specified."); - } - - - /* Do OMP! */ - - GAMMA_OUT = ompcoreGabor(D, x, DtX, 0, G, n, m, L, T, 0, gmode, profile, msgdelta, 0); - - return; -} -
--- a/util/ksvd utils/ompbox utils/ompmexGabor.m Thu Jul 21 16:37:14 2011 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,22 +0,0 @@ -%This is the Matlab interface to the OMP MEX implementation. -%The function is not for independent use, only through omp.m. - - -%OMPMEX Matlab interface to the OMP MEX implementation. -% GAMMA = OMPMEXGABOR(D,X,DtX,G,L,SPARSE_G,MSGDELTA,PROFILE) invokes the OMP -% MEX function according to the specified parameters. Not all the -% parameters are required. Those among D, X, DtX and G which are not -% specified should be passed as []. -% -% L - the target sparsity. -% SPARSE_G - returns a sparse GAMMA when nonzero, full GAMMA when zero. -% MSGDELTA - the delay in secs between messages. Zero means no messages. -% PROFILE - nonzero means that profiling information should be printed. - - -% Ron Rubinstein -% Computer Science Department -% Technion, Haifa 32000 Israel -% ronrubin@cs -% -% April 2009
--- a/util/ksvd utils/ompbox utils/ompprof.c Thu Jul 21 16:37:14 2011 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,113 +0,0 @@ -/************************************************************************** - * - * File name: ompprof.c - * - * Ron Rubinstein - * Computer Science Department - * Technion, Haifa 32000 Israel - * ronrubin@cs - * - * Last Updated: 11.4.2009 - * - *************************************************************************/ - - -#include "ompprof.h" - - -/* initialize profiling information */ - -void initprofdata(profdata *pd) -{ - pd->DtX_time = 0; - pd->XtX_time = 0; - pd->DtR_time = 0; - pd->maxabs_time = 0; - pd->DtD_time = 0; - pd->Lchol_time = 0; - pd->compcoef_time = 0; - pd->update_DtR_time = 0; - pd->update_resnorm_time = 0; - pd->compres_time = 0; - pd->indexsort_time = 0; - - pd->DtX_time_counted = 0; - pd->XtX_time_counted = 0; - pd->DtR_time_counted = 0; - pd->DtD_time_counted = 0; - pd->update_DtR_time_counted = 0; - pd->resnorm_time_counted = 0; - pd->compres_time_counted = 0; - pd->indexsort_time_counted = 0; - - pd->prevtime = clock(); -} - - -/* add elapsed time to profiling data according to specified computation */ - -void addproftime(profdata *pd, int comptype) -{ - switch(comptype) { - case DtX_TIME: pd->DtX_time += clock()-pd->prevtime; pd->DtX_time_counted = 1; break; - case XtX_TIME: pd->XtX_time += clock()-pd->prevtime; pd->XtX_time_counted = 1; break; - case DtR_TIME: pd->DtR_time += clock()-pd->prevtime; pd->DtR_time_counted = 1; break; - case DtD_TIME: pd->DtD_time += clock()-pd->prevtime; pd->DtD_time_counted = 1; break; - case COMPRES_TIME: pd->compres_time += clock()-pd->prevtime; pd->compres_time_counted = 1; break; - case UPDATE_DtR_TIME: pd->update_DtR_time += clock()-pd->prevtime; pd->update_DtR_time_counted = 1; break; - case UPDATE_RESNORM_TIME: pd->update_resnorm_time += clock()-pd->prevtime; pd->resnorm_time_counted = 1; break; - case INDEXSORT_TIME: pd->indexsort_time += clock()-pd->prevtime; pd->indexsort_time_counted = 1; break; - case MAXABS_TIME: pd->maxabs_time += clock()-pd->prevtime; break; - case LCHOL_TIME: pd->Lchol_time += clock()-pd->prevtime; break; - case COMPCOEF_TIME: pd->compcoef_time += clock()-pd->prevtime; break; - } - pd->prevtime = clock(); -} - - -/* print profiling info */ - -void printprofinfo(profdata *pd, int erroromp, int batchomp, int signum) -{ - clock_t tottime; - - tottime = pd->DtX_time + pd->XtX_time + pd->DtR_time + pd->DtD_time + pd->compres_time + pd->maxabs_time + - pd->Lchol_time + pd->compcoef_time + pd->update_DtR_time + pd->update_resnorm_time + pd->indexsort_time; - - mexPrintf("\n\n***** Profiling information for %s *****\n\n", erroromp? "OMP2" : "OMP"); - - mexPrintf("OMP mode: %s\n\n", batchomp? "Batch-OMP" : "OMP-Cholesky"); - - mexPrintf("Total signals processed: %d\n\n", signum); - - if (pd->DtX_time_counted) { - mexPrintf("Compute DtX time: %7.3lf seconds\n", pd->DtX_time/(double)CLOCKS_PER_SEC); - } - if (pd->XtX_time_counted) { - mexPrintf("Compute XtX time: %7.3lf seconds\n", pd->XtX_time/(double)CLOCKS_PER_SEC); - } - mexPrintf("Max abs time: %7.3lf seconds\n", pd->maxabs_time/(double)CLOCKS_PER_SEC); - if (pd->DtD_time_counted) { - mexPrintf("Compute DtD time: %7.3lf seconds\n", pd->DtD_time/(double)CLOCKS_PER_SEC); - } - mexPrintf("Lchol update time: %7.3lf seconds\n", pd->Lchol_time/(double)CLOCKS_PER_SEC); - mexPrintf("Compute coef time: %7.3lf seconds\n", pd->compcoef_time/(double)CLOCKS_PER_SEC); - if (pd->compres_time_counted) { - mexPrintf("Compute R time: %7.3lf seconds\n", pd->compres_time/(double)CLOCKS_PER_SEC); - } - if (pd->DtR_time_counted) { - mexPrintf("Compute DtR time: %7.3lf seconds\n", pd->DtR_time/(double)CLOCKS_PER_SEC); - } - if (pd->update_DtR_time_counted) { - mexPrintf("Update DtR time: %7.3lf seconds\n", pd->update_DtR_time/(double)CLOCKS_PER_SEC); - } - if (pd->resnorm_time_counted) { - mexPrintf("Update resnorm time: %7.3lf seconds\n", pd->update_resnorm_time/(double)CLOCKS_PER_SEC); - } - if (pd->indexsort_time_counted) { - mexPrintf("Index sort time: %7.3lf seconds\n", pd->indexsort_time/(double)CLOCKS_PER_SEC); - } - mexPrintf("---------------------------------------\n"); - mexPrintf("Total time: %7.3lf seconds\n\n", tottime/(double)CLOCKS_PER_SEC); -} -
--- a/util/ksvd utils/ompbox utils/ompprof.h Thu Jul 21 16:37:14 2011 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,106 +0,0 @@ -/************************************************************************** - * - * File name: ompprof.h - * - * Ron Rubinstein - * Computer Science Department - * Technion, Haifa 32000 Israel - * ronrubin@cs - * - * Last Updated: 18.8.2009 - * - * Collection of definitions and functions for profiling the OMP method. - * - *************************************************************************/ - - -#ifndef __OMP_PROF_H__ -#define __OMP_PROF_H__ - -#include "mex.h" -#include <time.h> - - - -/************************************************************************** - * - * Constants and data types. - * - **************************************************************************/ - - -/* constants denoting the various parts of the algorithm */ - -enum { DtX_TIME, XtX_TIME, DtR_TIME, MAXABS_TIME, DtD_TIME, LCHOL_TIME, COMPCOEF_TIME, - UPDATE_DtR_TIME, UPDATE_RESNORM_TIME, COMPRES_TIME, INDEXSORT_TIME }; - - - -/* profiling data container with counters for each part of the algorithm */ - -typedef struct profdata -{ - clock_t prevtime; /* the time when last initialization/call to addproftime() was performed */ - - clock_t DtX_time; - clock_t XtX_time; - clock_t DtR_time; - clock_t maxabs_time; - clock_t DtD_time; - clock_t Lchol_time; - clock_t compcoef_time; - clock_t update_DtR_time; - clock_t update_resnorm_time; - clock_t compres_time; - clock_t indexsort_time; - - /* flags indicating whether profiling data was gathered */ - int DtX_time_counted; - int XtX_time_counted; - int DtR_time_counted; - int DtD_time_counted; - int update_DtR_time_counted; - int resnorm_time_counted; - int compres_time_counted; - int indexsort_time_counted; - -} profdata; - - - -/************************************************************************** - * - * Initialize a profdata structure, zero all counters, and start its timer. - * - **************************************************************************/ -void initprofdata(profdata *pd); - - -/************************************************************************** - * - * Add elapsed time from last call to addproftime(), or from initialization - * of profdata, to the counter specified by comptype. comptype must be one - * of the constants in the enumeration above. - * - **************************************************************************/ -void addproftime(profdata *pd, int comptype); - - -/************************************************************************** - * - * Print the current contents of the counters in profdata. - * - * Parameters: - * pd - the profdata to print - * erroromp - indicates whether error-based (nonzero) or sparsity-based (zero) - * omp was performed. - * batchomp - indicates whether batch-omp (nonzero) or omp-cholesky (zero) - * omp was performed. - * signum - number of signals processed by omp - * - **************************************************************************/ -void printprofinfo(profdata *pd, int erroromp, int batchomp, int signum); - - -#endif -
--- a/util/ksvd utils/ompbox utils/omputils.c Thu Jul 21 16:37:14 2011 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,89 +0,0 @@ -/************************************************************************** - * - * File name: omputils.c - * - * Ron Rubinstein - * Computer Science Department - * Technion, Haifa 32000 Israel - * ronrubin@cs - * - * Last Updated: 18.8.2009 - * - *************************************************************************/ - -#include "omputils.h" -#include <math.h> - - -const char FULL_GAMMA_STR[] = "full"; -const char SPARSE_GAMMA_STR[] = "sparse"; - - -/* convert seconds to hours, minutes and seconds */ - -void secs2hms(double sectot, int *hrs, int *mins, double *secs) -{ - *hrs = (int)(floor(sectot/3600)+1e-2); - sectot = sectot - 3600*(*hrs); - *mins = (int)(floor(sectot/60)+1e-2); - *secs = sectot - 60*(*mins); -} - - -/* quicksort, public-domain C implementation by Darel Rex Finley. */ -/* modification: sorts the array data[] as well, according to the values in the array vals[] */ - -#define MAX_LEVELS 300 - -void quicksort(mwIndex vals[], double data[], mwIndex n) { - - long piv, beg[MAX_LEVELS], end[MAX_LEVELS], i=0, L, R, swap ; - double datapiv; - - beg[0]=0; - end[0]=n; - - while (i>=0) { - - L=beg[i]; - R=end[i]-1; - - if (L<R) { - - piv=vals[L]; - datapiv=data[L]; - - while (L<R) - { - while (vals[R]>=piv && L<R) - R--; - if (L<R) { - vals[L]=vals[R]; - data[L++]=data[R]; - } - - while (vals[L]<=piv && L<R) - L++; - if (L<R) { - vals[R]=vals[L]; - data[R--]=data[L]; - } - } - - vals[L]=piv; - data[L]=datapiv; - - beg[i+1]=L+1; - end[i+1]=end[i]; - end[i++]=L; - - if (end[i]-beg[i] > end[i-1]-beg[i-1]) { - swap=beg[i]; beg[i]=beg[i-1]; beg[i-1]=swap; - swap=end[i]; end[i]=end[i-1]; end[i-1]=swap; - } - } - else { - i--; - } - } -}
--- a/util/ksvd utils/ompbox utils/omputils.h Thu Jul 21 16:37:14 2011 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,77 +0,0 @@ -/************************************************************************** - * - * File name: omputils.h - * - * Ron Rubinstein - * Computer Science Department - * Technion, Haifa 32000 Israel - * ronrubin@cs - * - * Last Updated: 18.8.2009 - * - * Utility definitions and functions for the OMP library. - * - *************************************************************************/ - - -#ifndef __OMP_UTILS_H__ -#define __OMP_UTILS_H__ - -#include "mex.h" - - -/* constants for the representation mode of gamma */ - -extern const char FULL_GAMMA_STR[]; /* "full" */ -extern const char SPARSE_GAMMA_STR[]; /* "sparse" */ - - -#define FULL_GAMMA 1 -#define SPARSE_GAMMA 2 -#define INVALID_MODE 3 - - - -/************************************************************************** - * Memory management for OMP2. - * - * GAMMA_INC_FACTOR: - * The matrix GAMMA is allocated with sqrt(n)/2 coefficients per signal, - * for a total of nzmax = L*sqrt(n)/2 nonzeros. Whenever GAMMA needs to be - * increased, it is increased by a factor of GAMMA_INC_FACTOR. - * - * MAT_INC_FACTOR: - * The matrices Lchol, Gsub and Dsub are allocated with sqrt(n)/2 - * columns each. If additional columns are needed, this number is - * increased by a factor of MAT_INC_FACTOR. - **************************************************************************/ - -#define GAMMA_INC_FACTOR (1.4) -#define MAT_INC_FACTOR (1.6) - - - -/************************************************************************** - * Convert number of seconds to hour, minute and second representation. - * - * Parameters: - * sectot - total number of seconds - * hrs, mins, secs - output hours (whole) and minutes (whole) and seconds - * - **************************************************************************/ -void secs2hms(double sectot, int *hrs, int *mins, double *secs); - - - -/************************************************************************** - * QuickSort - public-domain C implementation by Darel Rex Finley. - * - * Modified to sort both the array vals[] and the array data[] according - * to the values in the array vals[]. - * - **************************************************************************/ -void quicksort(mwIndex vals[], double data[], mwIndex n); - - -#endif -
--- a/util/ksvd utils/ompbox utils/printf.m Thu Jul 21 16:37:14 2011 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,26 +0,0 @@ -function str = printf(varargin) -%PRINTF Print formatted text to screen. -% PRINTF(FMT,VAL1,VAL2,...) formats the data in VAL1,VAL2,... according to -% the format string FMT, and prints the result to the screen. -% -% The call to PRINTF(FMT,VAL1,VAL2,...) simply invokes the call -% DISP(SPRINTF(FMT,VAL1,VAL2,...)). For a complete description of the -% format string options see function SPRINTF. -% -% STR = PRINTF(...) also returns the formatted string. - - -% Ron Rubinstein -% Computer Science Department -% Technion, Haifa 32000 Israel -% ronrubin@cs -% -% April 2008 - - -if (nargout>0) - str = sprintf(varargin{:}); - disp(str); -else - disp(sprintf(varargin{:})); -end