plotcols( d,xcol,ycols,corr,colnames,LineSpecs)

function h = plotcols( d,xcol,ycols,corr,colnames,LineSpecs)
%PLOTCOLS plots data in columns of matrix d. One x column and any number of
% y columns can be chosen. Also the avg_on of the columns can automatically
% be subtracted and the columns can be divided by their standard deviation.
% This is very useful for comparing differently scaled data in one plot.
% Also, you can specify a name for each column, which is used for the
% legend and the x-axis label, and a LineSpec as in the plot function.
%
% If you find the help to long, try these commands:
%
%   plotcols(randn(1000,2)+100,2,0:2,'avgsd',{'C. F.' 'Gauss'},{'-' '.'})
%   plotcols(randn(1000,2)+100,2,0:2,'avg  ',{'C. F.' 'Gauss'},{'-' '.'})
%
%
% INPUTS:
% =======
%
% d:
% --
% d is a matrix with your x and y data in columns.
%
% xcol, ycols:
% ------------
% - x-axis data is taken from d(:,xcol). xcol is a scalar.
% - y-axis data is taken from d(:,ycols). ycols is a vector containing all
%   numbers of columns to be plot versus column xcol
% - an x or y column value of 0 corresponds to the row index, 1:N
%
% corr (optional):
% ----------------
% To get a comparable view of differently scaled data you can subtract
% average values of the column values, divide by the standard deviation or
% divide by the maximum absolute values (scale to -1:1).
% - set corr to 'avg' or 'avg' to subtract the average of each column
% - set corr to 'sd' to divide by the standard deviation
% - set corr to 'max' to divide columns by their max(abs) values
% - concatenate strings to do several operations (e.g. 'avgsd'). The
%   operation take place in the following order: avg -> sd -> max
%
% colnames (optional):
% --------------------
% To keep track of which columns are plotted, you can specify names for the
% columns in d. Those are displayed in the legend and in the x-axis of your
% plot
% - colnames is a a cell array of srings containing the name of each column
%   in d, e.g. colnames{3}='X' gives column 3 the label 'X'
% - the "0th column" is always named 'Row Index'
% - default names are 'column 1', 'column 2', ...
% 
% LineSpecs (optional):
% ---------------------
% You can specify a LineSpec for each column, just as in the plot command.
%
% - LineSpecs is a cell array of strings containing a LineSpec for each
%   column in d, or containing exactly one field, which is copied for all
%   columns. e.g. to plot column 2 dashed, set LineSpecs{2}='--', to plot
%   all columns dashed set LineSpecs={'--'}
% - default is a coloured solid line without a marker (same as in plot)
%
%
% OUTPUT:
% =======
% - h = plotcols(...) returns a column vector of handles to lineseries graphics
% objects, one handle per line.
%
%
% EXAMPLE:
% ========
%
% Suppose you have experimental data taken at certain points in time, e.g.
% time, temperature, x and y. This data is stored in a matrix:
%
%   d=[0.1 293.1 0.4 0.22
%      0.95 295.5 0.37 0.41
%      2.03 298.2 0.34 0.69
%      3.20 310.0 0.39 1.03
%      4.01 318.2 0.36 1.28
%      4.98 333.4 0.38 1.56]
%
% The column names are in a cell array:
%
%   c={'time' 'temperature' 'x' 'y'}
%
% Plot your data versus time with
%
%   plotcols(d,1,2:4,'',c)
%
% This does however not lead to a nice plot, since temperature has much
% higher numerical values than x and y. To see the plots on one level,
% subtract the avg_on value:
%
%   plotcols(d,1,2:4,'avg',c)
%
% Still the range of temperature is much higher than that of x and y.
% Divide by the standard deviation to fix that:
%
%   plotcols(d,1,2:4,'avgsd',c)
%
% Now you can see a nice correlation of T and y. Plot x and y versus T
% without avg and sd correction, also add markers:
%
%   plotcols(d,2,3:4,'',c,{'' '' 'rx-' 'bo-'})
%
% You might want to know how equally spaced time is. Plot time vs.
% row index:
%
%   plotcols(d,0,1,'',c)   or   plotcols(d,1,0,'',c)
%
%
% ----
% Version: 1.0, September 2008, Felix Zoergiebel (felix_z -> web.de)
%          - first version
% Version: 1.1, September 2008, Felix Zoergiebel (felix_z -> web.de)
%          - ignore both NaN and Inf when calculating averages
% Version: 1.2, November 2008, Felix Zoergiebel (felix_z -> web.de)
%          - added 'max' as an option for corr
%
% Feel free to use, modify and improve this script, but please always pass/
% update the version/author history above!
%
%
% See also plot, LineSpec

%% read corr 
if nargin<4 || isempty(corr), corr=''; end

avg_on=false; sd_on=false; max_on=false;
if ~isempty([findstr(corr,'mean') findstr(corr,'avg')]), avg_on=true; end
if ~isempty(findstr(corr,'sd')), sd_on=true; end
if ~isempty(findstr(corr,'max')), max_on=true; end

%% get dimensions, check ycols
ycols=ycols(:)';
ncolstot=size(d,2);
ncols=numel(ycols);
npoints=size(d,1);

if any(ycols<0 | ycols>ncolstot)
    error('All entries in ycols must be existing column numbers or 0.')
end

%% get and check column names and style
maxcol=max([xcol ycols]);
if nargin<5 || isempty(colnames)
    colnames=cell(maxcol,1);
    for idx=[ycols(ycols~=0) xcol(xcol~=0)], colnames{idx}=['column ' int2str(idx)]; end
elseif numel(colnames)<maxcol
    error('If you specify colnames, it must have at least an empty cell-entry for each column in xcol and ycols.')
else colnames=colnames(:);
end

if nargin<6, LineSpecs={''}; end
if ~iscell(LineSpecs), LineSpecs={LineSpecs}; end
if numel(LineSpecs)<=1
    LineSpecs=repmat(LineSpecs,max(ycols),1);
elseif numel(LineSpecs)<max(ycols)
    error('If you specify LineSpecs, it must have at least an empty cell-entry for each column in ycols OR exactly one cell-entry.')
else LineSpecs=LineSpecs(:);
end

%% get requested columns
if any([xcol ycols]==0)
    d=[(1:npoints)' d]; % set column 1 to 1:npoints
    xcol=xcol+1;
    ycols=ycols+1;
    colnames=[{'Row Index'};colnames]; % set 0th column's name to 'Row Index'
    LineSpecs=[{''};LineSpecs]; % default style for 0th column
end
d=d(:,[xcol ycols]); % extract requested columns

%% apply corrections
if avg_on || sd_on
    NaNInfidx=find(~isfinite(d(:,2:end))); % since it is just about a nice view just care about finite values.
    if avg_on
        dy0=d(:,2:end);
        dy0(NaNInfidx)=0;
        d(:,2:end)=bsxfun(@minus,d(:,2:end),sum(dy0,1)/npoints);
    end
    if sd_on
        dy0=d(:,2:end);
        dy0(NaNInfidx)=0;
        d(:,2:end)=bsxfun(@times,d(:,2:end),1./sqrt(sum(dy0.^2,1)/(npoints-1)));
    end
end
if max_on
    d(:,2:end)=bsxfun(@times,d(:,2:end),1./max(abs(d(:,2:end)),[],1));
end

%% setup plotting command
order=[ones(1,ncols);2:ncols+1];
plotcmd=mat2cell(d(:,order(:)),npoints,ones(1,2*ncols));
plotcmd=[reshape(plotcmd,2,ncols); LineSpecs(ycols)'];

%% plot data
if nargout>1
    h = plot(plotcmd{:});
else
    plot(plotcmd{:});
end

%% label graph if it goes to a new axis
if ~ishold
    
    title(['Columns from matrix ' inputname(1)])
    
    legend(colnames{ycols})
    
    xlabel(colnames{xcol})
    
    if ncols~=1, ylbl='original values';
    else ylbl=colnames{ycols};
    end
    if avg_on, ylbl=[ylbl ' - AVGs'];
        if sd_on || max_on, ylbl=['( ' ylbl ' )']; end
    end
    if sd_on, ylbl=[ylbl ' / SDs']; end
    if max_on, ylbl=[ylbl ' / max. abs.']; end
    ylabel(ylbl)
    
end

end