function hL = dsplot(x, y, numPoints)

%DSPLOT Create down sampled plot.
%   This function creates a down sampled plot to improve the speed of
%   exploration (zoom, pan).
%
%   DSPLOT(X, Y) plots Y versus X by downsampling if there are large number
%   of elements. X and Y needs to obey the following:
%     1. X must be a monotonically increasing vector.
%     2. If Y is a vector, it must be the same size as X.
%     3. If Y is a matrix, one of the dimensions must line up with X.
%
%   DSPLOT(Y) plots the columns of Y versus their index.
%
%   hLine = DSPLOT(X, Y) returns the handles of the line. Note that the
%   lines may be downsampled, so they may not represent the full data set.
%
%   DSPLOT(X, Y, NUMPOINTS) or DSPLOT(Y, [], NUMPOINTS) specifies the
%   number of points (roughly) to display on the screen. The default is
%   50000 points (~390 kB doubles). NUMPOINTS can be a number greater than
%   500.
%
%   It is very likely that more points will be displayed than specified by
%   NUMPOINTS, because it will try to plot any outlier points in the range.
%   If the signal is stochastic or has a lot of sharp changes, there will
%   be more points on plotted on the screen.
%
%   The figure title (name) will indicate whether the plot shown is
%   downsampled or is the true representation.
%
%   The figure can be saved as a .fig file, which will include the actual
%   data. The figure can be reloaded and the actual data can be exported to
%   the base workspace via a menu.
%
%   Run the following examples and zoom/pan to see the performance.
%
%  Example 1: (with small details)
%   x  = linspace(0, 2*pi, 1000000);
%   y1 = sin(x)+.02*cos(200*x)+0.001*sin(2000*x)+0.0001*cos(20000*x);
%   dsplot(x,y1);title('Down Sampled');
%   % compare with
%   figure;plot(x,y1);title('Normal Plot');
%
%  Example 2: (with outlier points)
%   x  = linspace(0, 2*pi, 1000000);
%   y1 = sin(x) + .01*cos(200*x) + 0.001*sin(2000*x);
%   y2 = sin(x) + 0.3*cos(3*x)   + 0.001*randn(size(x));
%   y1([300000, 700000, 700001, 900000]) = [0, 1, -2, 0.5];
%   y2(300000:500000) = y2(300000:500000) + 1;
%   y2(500001:600000) = y2(500001:600000) - 1;
%   y2(800000) = 0;
%   dsplot(x, [y1;y2]);title('Down Sampled');
%   % compare with
%   figure;plot(x, [y1;y2]);title('Normal Plot');
%
%  See also PLOT.

%  Version:
%   v1.0 - first version (Aug 1, 2007)
%   v1.1 - added CreateFcn for the figure so that when the figure is saved
%          and re-loaded, the zooming and panning works. Also added a menu
%          item for saving out the original data back to the base
%          workspace. (Aug 10, 2007)
%
%  Jiro Doke
%  August 1, 2007

debugMode = false;

%--------------------------------------------------------------------------
% Error checking
error(nargchk(1, 3, nargin, 'struct'));
if nargin < 3
  % Number of points to show on the screen. It's quite possible that more
  % points will be displayed if there are outlier points
  numPoints = 50000;  % ~390 kB for doubles
end
if nargin == 1 || isempty(y)
  noXVar = true;
  y = x;
  x = [];
else
  noXVar = false;
end
myErrorCheck;
%--------------------------------------------------------------------------

if size(x, 2) > 1  % it's a row vector -> transpose
  x = x';
  y = y';
  varTranspose = true;
else
  varTranspose = false;
end

% Number of lines
numSignals = size(y, 2);

% If the number of lines is greater than the number of data points per
% line, it's possible that the user may have mistaken the matrix
% orientation.
if numSignals > size(y, 1)
  s = input(sprintf('Are you sure you want to plot %d lines? (y/n) ', ...
    numSignals), 's');
  if ~strcmpi(s, 'y')
    disp('Canceled. You may want to transpose the matrix.');
    if nargout == 1
      hL = [];
    end
    return;
  end
end

% Attempt to find outliers. Use a running average technique
filterWidth = ceil(min([50, length(x)/10])); % max window size of 50
a  = y - filter(ones(filterWidth,1)/filterWidth, 1, y);
[iOutliers, jOutliers] = find(abs(a - repmat(mean(a), size(a, 1), 1)) > ...
  repmat(4 * std(a), size(a, 1), 1));
clear a;

% Always create new figure because it messes around with zoom, pan,
% datacursors.
hFig    = figure;
figName = '';

% Create template plot using NaNs
hLine   = plot(NaN(2, numSignals), NaN(2, numSignals));
set(hLine, 'tag', 'dsplot_lines');

% Define CreateFcn for the figure
set(hFig, 'CreateFcn', @mycreatefcn);
mycreatefcn();

% Create menu for exporting data
hMenu = uimenu(hFig, 'Label', 'Data');
uimenu(hMenu, ...
  'Label'   , 'Export data to workspace.', ...
  'Callback', @myExportFcn);

% Update lines
updateLines([min(x), max(x)]);

% Deal with output argument
if nargout == 1
  hL = hLine;
end

%--------------------------------------------------------------------------
  function myExportFcn(varargin)
    % This callback allows for extracting the actual data from the figure.
    % This means that if you save this figure and load it back later, you
    % can get back the data.
    
    % Determine the variable name
    allVarNames = evalin('base', 'who');
    newVarName = genvarname('dsplotData', allVarNames);
    
    % X
    if ~noXVar
      if varTranspose
        dat.x = x';
      else
        dat.x = x;
      end
    end
    
    % Y
    if varTranspose
      dat.y = y';
    else
      dat.y = y;
    end
    
    assignin('base', newVarName, dat);
    
    msgbox(sprintf('Data saved to the base workspace as ''%s''.', ...
      newVarName), 'Saved', 'modal');
    
  end

%--------------------------------------------------------------------------
  function mycreatefcn(varargin)
    % This callback defines the custom zoom/pan functions. It is defined as
    % the CreateFcn of the figure, so it allows for saving and reloading of
    % the figure.

    if nargin > 0
      hFig = varargin{1};
    end
    hLine = findobj(hFig, 'type', 'axes');
    hLine(strmatch('legend', get(hLine, 'tag'))) = [];
    hLine = get(hLine, 'Children');
    
    % Create Zoom, Pan, Datacursor objects
    hZoom = zoom(hFig);
    hPan  = pan(hFig);
    hDc   = datacursormode(hFig);
    set(hZoom, 'ActionPostCallback', @mypostcallback);
    set(hPan , 'ActionPostCallback', @mypostcallback);
    set(hDc  , 'UpdateFcn'         , @myDCupdatefcn);

  end

%--------------------------------------------------------------------------
  function mypostcallback(obj, evd) %#ok
    % This callback that gets called when the mouse is released after
    % zooming or panning.

    % single or double-click
    switch get(hFig, 'SelectionType')
      case {'normal', 'alt'}
        updateLines(xlim(evd.Axes));

      case 'open'
        updateLines([min(x), max(x)]);

    end

  end

%--------------------------------------------------------------------------
  function updateLines(rng)
    % This helper function is for determining the points to plot on the
    % screen based on which portion is visible in the current limits.

    % find indeces inside the range
    id = find(x >= rng(1) & x <= rng(2));

    % if there are more points than we want
    if length(id) > numPoints / numSignals

      % see how many outlier points are in this range
      blah = iOutliers > id(1) & iOutliers < id(end);

      % determine indeces of points to plot. 
      idid = round(linspace(id(1), id(end), round(numPoints/numSignals)))';

      x2 = cell(numSignals, 1);
      y2 = x2;
      for iSignals = 1:numSignals
        % add outlier points
        ididid = unique([idid; iOutliers(blah & jOutliers == iSignals)]);
        x2{iSignals} = x(ididid);
        y2{iSignals} = y(ididid, iSignals);
      end

      if debugMode
        figName = ['downsampled - ', sprintf('%d, ', cellfun('length', y2))];
      else
        figName = 'downsampled';
      end

    else % no need to down sample
      figName = 'true';

      x2 = repmat({x(id)}, numSignals, 1);
      y2 = mat2cell(y(id, :), length(id), ones(1, numSignals))';

    end

    % Update plot
    set(hLine, {'xdata', 'ydata'} , [x2, y2]);
    set(hFig, 'Name', figName);

  end

%--------------------------------------------------------------------------
  function txt = myDCupdatefcn(empt, event_obj) %#ok
    % This function displays appropriate data cursor message based on the
    % display type

    pos = get(event_obj,'Position');
    switch figName
      case 'true'
        txt = {['X: ',num2str(pos(1))],...
          ['Y: ',num2str(pos(2))]};
      otherwise
        txt = {['X: ',num2str(pos(1))],...
          ['Y: ',num2str(pos(2))], ...
          'Warning: Downsampled', ...
          'May not be accurate'};
    end
  end

%--------------------------------------------------------------------------
  function myErrorCheck
    % Do some error checking on the input arguments.

    if ~isa(numPoints, 'double') || numel(numPoints) > 1 || numPoints < 500
      error('Third argument must be a scalar greater than 500');
    end
    if ~isnumeric(x) || ~isnumeric(y)
      error('Arguments must be numeric');
    end
    if length(size(x)) > 2 || length(size(y)) > 2
      error('Only 2-D data accepted');
    end
    
    % If only one input, create index vector X
    if isempty(x)
      if ismember(1, size(y))
        x = reshape(1:numel(y), size(y));
      else
        x = (1:size(y, 1))';
      end
    end
    
    if ~ismember(1, size(x))
      error('First argument has to be a vector');
    end
    if ~isequal(size(x, 1), size(y, 1)) && ~isequal(size(x, 2), size(y, 2))
      error('One of the dimensions of the two arguments must match');
    end
    if any(diff(x) <= 0)
      error('The first argument has to be a monotonically increasing vector');
    end
  end

end