ivan@140: function gamma = omp2Gabor(varargin)
ivan@140: %% omp2Gabor Error-constrained Orthogonal Matching Pursuit.
ivan@140: %
ivan@140: %  This file is part of SMALLbox [2] and it is adaptation of Ron
ivan@140: %  Rubenstein omp solver [1] for Gabor dictionary as defined in 
ivan@140: %  Audio Inpainting by Adler et al [3]. The dictionary is presented as
ivan@140: %  DCT+DST and solver pics two atoms per iteration (given the frequency one
ivan@140: %  from DCT and one from DST). For more information look in [3].
ivan@140: %  
ivan@140: %  GAMMA = OMP2(D,X,G,EPSILON) solves the optimization problem
ivan@140: %
ivan@140: %       min  |GAMMA|_0     s.t.  |X - D*GAMMA|_2 <= EPSILON
ivan@140: %      gamma
ivan@140: %
ivan@140: %  for each of the signals in X, using Batch Orthogonal Matching Pursuit.
ivan@140: %  Here, D is a dictionary with normalized columns, X is a matrix
ivan@140: %  containing column signals, EPSILON is the error target for each signal,
ivan@140: %  and G is the Gramm matrix D'*D. The output GAMMA is a matrix containing
ivan@140: %  the sparse representations as its columns. 
ivan@140: %
ivan@140: %  GAMMA = OMP2(D,X,[],EPSILON) performs the same operation, but without
ivan@140: %  the matrix G, using OMP-Cholesky. This call produces the same output as
ivan@140: %  Batch-OMP, but is significantly slower. Using this syntax is only
ivan@140: %  recommended when available memory is too small to store G.
ivan@140: %
ivan@140: %  GAMMA = OMP2(DtX,XtX,G,EPSILON) is the fastest implementation of OMP2,
ivan@140: %  but also requires the most memory. Here, DtX stores the projections
ivan@140: %  D'*X, and XtX is a row vector containing the squared norms of the
ivan@140: %  signals, sum(X.*X). In this case Batch-OMP is used, but without having
ivan@140: %  to compute D'*X and XtX in advance, which slightly improves runtime.
ivan@140: %  Note that in general, the call
ivan@140: %
ivan@140: %    GAMMA = OMP2(D'*X, sum(X.*X), G, EPSILON);
ivan@140: %
ivan@140: %  will be faster than the call
ivan@140: %
ivan@140: %    GAMMA = OMP2(D,X,G,EPSILON);
ivan@140: %
ivan@140: %  due to optimized matrix multiplications in Matlab. However, when the
ivan@140: %  entire matrix D'*X cannot be stored in memory, one of the other two
ivan@140: %  versions can be used. Both compute D'*X for just one signal at a time,
ivan@140: %  and thus require much less memory.
ivan@140: %
ivan@140: %  GAMMA = OMP2(...,PARAM1,VAL1,PARAM2,VAL2,...) specifies additional
ivan@140: %  parameters for OMP2. Available parameters are:
ivan@140: %
ivan@140: %    'gammamode' - Specifies the representation mode for GAMMA. Can be
ivan@140: %                  either 'full' or 'sparse', corresponding to a full or
ivan@140: %                  sparse matrix, respectively. By default, GAMMA is
ivan@140: %                  returned as a sparse matrix.
ivan@140: %    'maxatoms' -  Limits the number of atoms in the representation of each
ivan@140: %                  signal. If specified, the number of atoms in each
ivan@140: %                  representation does not exceed this number, even if the
ivan@140: %                  error target is not met. Specifying maxatoms<0 implies
ivan@140: %                  no limit (default).
ivan@140: %    'messages'  - Specifies whether progress messages should be displayed.
ivan@140: %                  When positive, this is the number of seconds between
ivan@140: %                  status prints. When negative, indicates that no messages
ivan@140: %                  should be displayed (this is the default).
ivan@140: %    'checkdict' - Specifies whether dictionary normalization should be
ivan@140: %                  verified. When set to 'on' (default) the dictionary
ivan@140: %                  atoms are verified to be of unit L2-norm. Setting this
ivan@140: %                  parameter to 'off' disables verification and accelerates
ivan@140: %                  function performance. Note that an unnormalized
ivan@140: %                  dictionary will produce invalid results.
ivan@140: %    'profile'   - Can be either 'on' or 'off'. When 'on', profiling
ivan@140: %                  information is displayed at the end of the funciton
ivan@140: %                  execution.
ivan@140: %
ivan@140: %
ivan@140: %  Summary of OMP2 versions:
ivan@140: %
ivan@140: %    version                 |   speed     |   memory
ivan@140: %  -------------------------------------------------------------
ivan@140: %   OMP2(DtX,XtX,G,EPSILON)  |  very fast  |  very large
ivan@140: %   OMP2(D,X,G,EPSILON)      |  fast       |  moderate
ivan@140: %   OMP2(D,X,[],EPSILON)     |  very slow  |  small
ivan@140: %  -------------------------------------------------------------
ivan@140: %
ivan@140: %
ivan@140: %  References:
ivan@140: %  [1] M. Elad, R. Rubinstein, and M. Zibulevsky, "Efficient Implementation
ivan@140: %      of the K-SVD Algorithm using Batch Orthogonal Matching Pursuit",
ivan@140: %      Technical Report - CS, Technion, April 2008.
ivan@140: %
ivan@140: %  See also OMP.
ivan@140: 
ivan@140: 
ivan@140: %  Ron Rubinstein
ivan@140: %  Computer Science Department
ivan@140: %  Technion, Haifa 32000 Israel
ivan@140: %  ronrubin@cs
ivan@140: %
ivan@140: %  April 2009
ivan@140: 
ivan@140: %
ivan@140: %   Centre for Digital Music, Queen Mary, University of London.
ivan@140: %   This file copyright 2011 Ivan Damnjanovic.
ivan@140: %
ivan@140: %   This program is free software; you can redistribute it and/or
ivan@140: %   modify it under the terms of the GNU General Public License as
ivan@140: %   published by the Free Software Foundation; either version 2 of the
ivan@140: %   License, or (at your option) any later version.  See the file
ivan@140: %   COPYING included with this distribution for more information.
ivan@140: %%
ivan@140: 
ivan@140: % default options
ivan@140: 
ivan@140: sparse_gamma = 1;
ivan@140: msgdelta = -1;
ivan@140: maxatoms = -1;
ivan@140: checkdict = 1;
ivan@140: profile = 0;
ivan@140: 
ivan@140: 
ivan@140: % determine number of parameters
ivan@140: 
ivan@140: paramnum = 1;
ivan@140: while (paramnum<=nargin && ~ischar(varargin{paramnum}))
ivan@140:   paramnum = paramnum+1;
ivan@140: end
ivan@140: paramnum = paramnum-1;
ivan@140: 
ivan@140: 
ivan@140: % parse options
ivan@140: 
ivan@140: for i = paramnum+1:2:length(varargin)
ivan@140:   paramname = varargin{i};
ivan@140:   paramval = varargin{i+1};
ivan@140: 
ivan@140:   switch lower(paramname)
ivan@140: 
ivan@140:     case 'gammamode'
ivan@140:       if (strcmpi(paramval,'sparse'))
ivan@140:         sparse_gamma = 1;
ivan@140:       elseif (strcmpi(paramval,'full'))
ivan@140:         sparse_gamma = 0;
ivan@140:       else
ivan@140:         error('Invalid GAMMA mode');
ivan@140:       end
ivan@140:       
ivan@140:     case 'maxatoms'
ivan@140:       maxatoms = paramval;
ivan@140:       
ivan@140:     case 'messages'
ivan@140:       msgdelta = paramval;
ivan@140: 
ivan@140:     case 'checkdict'
ivan@140:       if (strcmpi(paramval,'on'))
ivan@140:         checkdict = 1;
ivan@140:       elseif (strcmpi(paramval,'off'))
ivan@140:         checkdict = 0;
ivan@140:       else
ivan@140:         error('Invalid checkdict option');
ivan@140:       end
ivan@140: 
ivan@140:     case 'profile'
ivan@140:       if (strcmpi(paramval,'on'))
ivan@140:         profile = 1;
ivan@140:       elseif (strcmpi(paramval,'off'))
ivan@140:         profile = 0;
ivan@140:       else
ivan@140:         error('Invalid profile mode');
ivan@140:       end
ivan@140: 
ivan@140:     otherwise
ivan@140:       error(['Unknown option: ' paramname]);
ivan@140:   end
ivan@140:   
ivan@140: end
ivan@140: 
ivan@140: 
ivan@140: % determine call type
ivan@140: 
ivan@140: if (paramnum==4)
ivan@140:   
ivan@140:   n1 = size(varargin{1},1);
ivan@140:   n2 = size(varargin{2},1);
ivan@140:   n3 = size(varargin{3},1);
ivan@140:   
ivan@140:   if ( (n1>1 && n2==1) || (n1==1 && n2==1 && n3==1) )  %  DtX,XtX,G,EPSILON
ivan@140:     
ivan@140:     DtX = varargin{1};
ivan@140:     XtX = varargin{2};
ivan@140:     G = varargin{3};
ivan@140:     epsilon = varargin{4};
ivan@140:     D = [];
ivan@140:     X = [];
ivan@140:     
ivan@140:   else  % D,X,G,EPSILON
ivan@140:     
ivan@140:     D = varargin{1};
ivan@140:     X = varargin{2};
ivan@140:     G = varargin{3};
ivan@140:     epsilon = varargin{4};
ivan@140:     DtX = [];
ivan@140:     XtX = [];
ivan@140:     
ivan@140:   end
ivan@140:   
ivan@140: else
ivan@140:   error('Invalid number of parameters');
ivan@140: end
ivan@140: 
ivan@140: G=[];
ivan@140: 
ivan@140: % verify dictionary normalization
ivan@140: 
ivan@140: if (checkdict)
ivan@140:   if (isempty(G))
ivan@140:     atomnorms = sum(D.*D);
ivan@140:   else
ivan@140:     atomnorms = diag(G);
ivan@140:   end
ivan@140:   if (any(abs(atomnorms-1) > 1e-2))
ivan@140:     error('Dictionary columns must be normalized to unit length');
ivan@140:   end
ivan@140: end
ivan@140: 
ivan@140: 
ivan@140: % omp
ivan@140: 
ivan@140: gamma = omp2mexGabor(D,X,DtX,XtX,G,epsilon,sparse_gamma,msgdelta,maxatoms,profile);