function Y = mslowess(MZ,Y,varargin)
%MSLOWESS provides nonparametric smoothing of mass spectrometry signals
%
%   YOUT = MSLOWESS(MZ,Y) smoothes the mass spectrometry signal Y using the
%   'Lowess' method with a default span of 10 samples. MZ is the
%   mass-charge column vector with the same size as the vector of ion
%   intensity samples Y (or spectrum). Y can be a matrix with several
%   spectrograms, all sharing the same MZ scale.
%
%   Notes:
%   1) MSLOWESS assumes an MZ vctor that is not uniformly spaced; therefore,
%   the sliding window for smoothing is centered using the closest samples
%   in terms of the MZ value and not in terms of the MZ index.
%   2) When the MZ vector does not have repeated values or NaNs, the
%   algorithm will be approximately two times faster.
%
%   MSLOWESS(...,'ORDER',O) sets the order O of the 'lowess' smoother.
%   Order can be 1 (linear fit or Lowess) or 2 (quadratic fit or Loess).
%   The default order is 1. Order can also be 0, which is equivalent to a
%   weighted local mean estimator, presumably faster, because only a mean
%   computation is performed instead of a least squares regression. 
%
%   Note: the MATLAB Curve Fitting Toolbox refers to Lowess smoothing of
%   order 2 as Loess smoothing. 
%
%   MSLOWESS(...,'SPAN',S) modifies the window size for the smoothing
%   kernel. If SPAN is greater than 1, then the window is of size span (in
%   samples, i.e., independently of the MZ vector). The default value is 10
%   samples. Higher values smooth the signal more at the expense of
%   computation time.  If S is less than 1, then the window size is taken
%   to be a fraction of the number of points in the data; e.g., for S =
%   0.005, the window size is equal to 0.50% of the number of points in MZ.
%
%   MSLOWESS(...,'KERNEL',K) selects the kernel function. The kernel
%   function is used to weight the observed ion intensities such that
%   those samples close to the MZ location being smoothed have most weight
%   in determining the estimate. Options are
%
%       'tricubic' (default)    (1 - (dist/dmax).^3).^3  
%       'gaussian'              exp(-(2*dist/dmax).^2)
%       'linear'                1-dist/dmax
%
%   MSLOWESS(...,'ROBUSTITERATIONS',I) sets the number of iterations of a
%   robust fit. If I is 0 (default), no robust fit is performed. For robust
%   smoothing, small residual values at every span are weighted to
%   improve the new estimate. 1 or 2 robust iterations are usually
%   adequate; larger values can be computationally expensive.
%
%   Note: For a uniformly spaced MZ vector a nonrobust smoothing with
%   order 0 is equivalent to filtering the signal with the KERNEL vector.
%
%   MSLOWESS(...,'SHOWPLOT',SP) plots the smoothed spectra over the
%   original. When SP is TRUE, the first spectrogram in Y is used. If
%   MSLOWESS is called without output arguments, a plot will be shown
%   unless SP is FALSE. SP can also contain an index to one of the
%   spectrograms in Y.
%
%   Example: 
%
%       load sample_lo_res
%
%       % Smooth a group of spectrograms.
%       YS = mslowess(MZ_lo_res,Y_lo_res);
%
%       % Plot the third spectrogram in Y_lo_res and its smoothed signal.
%       mslowess(MZ_lo_res,Y_lo_res,'SHOWPLOT',3);
%
%   See also MSALIGN, MSBACKADJ, MSHEATMAP, MSNORM, MSRESAMPLE, MSSGOLAY,
%   MSVIEWER. 

% Copyright 2003-2004 The MathWorks, Inc.
% $Revision: 1.1.12.1 $   $Date: 2004/12/24 20:43:27 $


% set default parameters
kernel = 'tricubic';
span = 10;
order = 1; % 'lowess'
robustIter = 0;
if nargout == 1
    plotId = 0; 
else
    plotId = 1;
end


% deal with the various inputs
if nargin > 2
    if rem(nargin,2) == 1
        error('Bioinfo:mslowess:IncorrectNumberOfArguments',...
            'Incorrect number of arguments to %s.',mfilename);
    end
    okargs = {'kernel','robust','span','order','showplot'};
    for j=1:2:nargin-2
        pname = varargin{j};
        pval = varargin{j+1};
        k = strmatch(lower(pname), okargs);%#ok
        if isempty(k)
            error('Bioinfo:mslowess:UnknownParameterName',...
                'Unknown parameter name: %s.',pname);
        elseif length(k)>1
            error('Bioinfo:mslowess:AmbiguousParameterName',...
                'Ambiguous parameter name: %s.',pname);
        else
            switch(k)
                case 1  % method
                    kernels = {'tricubic','gaussian','linear'};
                    kernel = strmatch(lower(pval),kernels); %#ok
                    if isempty(kernel)
                        error('Bioinfo:mslowess:NotValidKernel',...
                              'Not a valid smoothing kernel.')
                    end
                    kernel=kernels{kernel};
                case 2 % robust iterations
                    if ~isnumeric(pval) || ~isscalar(pval)
                        error('Bioinfo:mslowess:LowessSpanNumeric',...
                              'SPAN must be a numeric value.');
                    end
                    robustIter = pval;
                    if robustIter < 0
                        error('Bioinfo:mslowess:IterationsPositive',...
                              'ITERATIONS must be a positive value.');
                    end
                case 3  % span
                    if ~isnumeric(pval) || ~isscalar(pval)
                        error('Bioinfo:mslowess:LowessSpanNumeric',...
                            'SPAN must be a numeric value.');
                    end
                    span = pval;
                    if span < 0
                        error('Bioinfo:mslowess:LowessSpanPositive',...
                            'SPAN must be a positive value.');
                    end
                    if span >numel(MZ)
                        error('Bioinfo:mslowess:LowessSpanLessThanX',...
                              'SPAN must be less than the number of elements in Y.');
                    end
                case 4 % order
                    if ~isnumeric(pval) || ~isscalar(pval)
                        error('Bioinfo:mslowess:OrderNumeric',...
                              'ORDER must be a numeric value.');
                    end
                    order = pval;
                    if ~ismember(order,[0 1 2])
                        error('Bioinfo:mslowess:LowessBadNumericOrder',...
                              'ORDER should be 0, 1, or 2.');
                    end
                 case 5 % show
                    if opttf(pval) 
                        if isnumeric(pval)
                            if isscalar(pval)
                                plotId = double(pval); 
                            else
                                plotId = 1;
                                warning('Bioinfo:mslowess:SPNoScalar',...
                                'SHOWPLOT must be a single index to one of the spectrograms in Y or\na logical. Plotting the first column in Y only.')
                            end
                        else
                            plotId = 1;
                        end
                    else
                        plotId = 0;
                    end
            end
        end
    end
end

% set the method to call masmooth
methods = {'mean','lowess','loess'};
method = methods{order + 1};
if robustIter > 0
    method = ['r',method];
end

% validate MZ and Y
if size(MZ,2) ~= 1
   error('Bioinfo:mslowess:OnlyOneMZScale','MZ scale must be a single column.')
end
if ~isnumeric(Y) || ~isreal(Y)
   error('Bioinfo:mslowess:IntensityNotNumericAndReal',...
         'Ion intensities must be numeric and real.') 
end
if ~isnumeric(MZ) || ~isreal(MZ)
   error('Bioinfo:mslowess:MZNotNumericAndReal','MZ scale must be numeric and real.') 
end
if size(MZ,1) ~= size(Y,1)
   error('Bioinfo:mslowess:NotEqualNumberOfSamples',...
         ['The ion intensity vector(s) must have the\n'...
          'same number of samples as the MZ scale.'])
end

numSpectrograms = size(Y,2);

if (plotId~=0) && ~any((1:numSpectrograms)==plotId)
    warning('Bioinfo:mslowess:InvalidPlotIndex',...
    'SHOWPLOT is not a valid spectrogram index, no plot was generated.')
end

% iterate for every spectrogram
for ns = 1:numSpectrograms
if nargout>0 || (ns == plotId)
    
    if (ns == plotId)
        yo = Y(:,ns);
    end
    
    Y(:,ns) = masmooth(MZ,Y(:,ns),span,method,robustIter,kernel);
    
    if (ns == plotId)
        figure
        plot(MZ,[yo Y(:,ns)])
        title(sprintf('Spectrogram ID: %d',ns));
        xlabel('Mass/Charge (M/Z)')
        ylabel('Relative Intensity')
        legend('Original spectrogram','Smoothed spectrogram')
        axis([min(MZ) max(MZ) min(Y(:,ns)) max(Y(:,ns))])
        grid on
        hold off
    end
    
end % if nargout>0 || (ns == plotId)    
end % for ns = 1:numSpectrograms