meansemfun(datasource, xcolindex, ycolindex, paramcolindex)

function [curvestructarray, relevant_data]  = meansemfun(datasource, xcolindex, ycolindex, paramcolindex)
% function curvestructarray =  meansemfun(datasource, xcolindex, ycolindex,
% paramcolindex)
% Meansemfun plots y vs x with errorbars, starting with a disordered list
% of data. The data should be in the form of a matrix where each row
% represents an observation, with columns for the dependent variable (y)
% and for an independent (x axis) variable. Optionally it can also have a
% column for a second independent variable, that is used for distinguishing
% among different curves to be plotted. There can be any number of rows,
% with x and curve-parameter values recurring in any order and with any
% frequency. The x values can be different for the different curves.
%
% If no arguments are supplied, the user is asked to specify a file source
% for the data, by entering the name of a text file that contains the data
% matrix. Arguments 2:4 are the column indices for x, y, and (optionally)
% the curve parameter; if these are missing, user input is requested.
% Meansemfun plots the mean and +-1 x the s.e.m. of y as a function of x
% for each value of the curve parameter and returns curvestructarray, an
% array of structures, one for each curve. Fields are X; mean of Y;
% standard error of the mean of Y; and number of observations, each as a
% column vectors ordered in ascending order of X; and the curve parameter
% value, if applicable. Also returned are the relevant data, suitably
% sorted.
% EXAMPLE:
%
% If datalist is a matrix in the workspace (primary usage):
% >>meansemfun(datalist,1,2) creates an errorbar plot.
% Y values are taken from datalist column 2
% X values are taken from datalist column 1
%
% >>meansemfun(datalist,1,2,4) creates a plot containing one errorbar curve
% for each of the Z (column 4) values that appear in column 4 of datalist.
% Each curve is based on all the data for the corresponding Z value.
% Y axis shows mean and standard error of Y values from datalist column 2;
% X axis shows the corresponding X values from datalist column 1;
% legend shows value(s) of Z, (datalist column 4) that apply to each curve.
%
% To load the data from an ascii file in the current directory:
% >>meansemfun('mydatalist.dat') loads a data matrix from that file.
% You are then prompted for the column indices for x, y, and optionally,
% if there are many curves to plot, z. Or, you can use
% >>meansemfun('mydatalist.dat', 1,2,5)
% if you already know that x is column 1, and y and z are columns 2 and 5.
% >>meansemfun
% will prompt the user for an ascii data file name and column indices.
%
% mycurvestructs = meansemfun... returns an array of structures, one for
% each value of Z. Each structure has these fields:
%     x     % vector of all x values that occur for this z value 
%     ymean  % mean of y for corresponding z and x
%     ysem
%     nmeasures % the number of replications (data rows) for given z and x.
%     z  % (the additional independent variable value for that curve)
%     plotsym  % the plot symbol, color and line style for that curve
% Here mycurvestructs(1) has information for the lowest-Z curve, etc.
% Plot symbols cycle through 4 symbols and colors, with continuous lines,
% and then the same sequence with dashed lines, repeating if > 8 curves.

% don macleod, ed vul, 1993, 2008 
% July 2008: less clunky code, color graphics

if nargin==2 || nargin>4
    error('meansemfun(datasource, xcolindex, ycolindex, paramcolindex) requires 0,1,3 or 4 arguments');
end
%% Ask user for data source and relevant column indices as needed:
if nargin==0 % have to ask for data file as source
% File assumed to be in current directory; consider uigetfile()
    datasourcename = input('ascii file to analyze, including any extension: ','s');
    load(datasourcename); % creates a data matrix variable, named with the filename prefix
    if any(datasourcename==46)  % found an extension (46 = ascii('.'))
        datasourcename = datasourcename(1:find(datasourcename==46)-1); % trim off the extension
    end
    originaldata = eval(datasourcename); % datasourcename string has name of loaded data matrix variable
else   % data matrix is 1st argument
    originaldata = datasource;
end
if isempty(originaldata)
    error('No data in specified source');
end
if nargin<4
    paramcolindex = [];
end
if nargin<=1
    xcolindex = input('Column index (>=1) for x: ');
    ycolindex = input('Column index for y: ');
    paramcolindex = input(...
        'Column index for z, to create multiple curves, one for each z value (hit ''enter'' if none)? ');
end
if isempty(paramcolindex)                       % want just a single curve
    relevant_data = [originaldata(:,[xcolindex ycolindex]) ones(size(originaldata,1),1) ];  % [X Y Z] as 3 columns, Z is dummy
else
    relevant_data = originaldata(:,[xcolindex ycolindex paramcolindex]);
end
%% sort data, get means and standard errors for y per z value per x value
zs = sort(unique(relevant_data(:,3)));
xs = cell(size(zs));
figure();
% cell array containing strings for plot colors and symbols:
plotsym = {'-k*' '-rd' '-go' '-bs'  '-m+' ':rd' '-yo' ':gs' ':bo'   };  
for zi = 1:length(zs)
    legs{zi} = ['Z = ' num2str(zs(zi))];
    xydata = relevant_data((relevant_data(:,3)==zs(zi)),1:2); 
    xs{zi}  = sort(unique(xydata(:,1)));
    ymean = [];
    ysem = [];
    for xi = 1:length(xs{zi})
        ydata = xydata((xydata(:,1)==xs{zi}(xi)),2);
        ymean(xi) = mean(ydata);
        ysem(xi) = std(ydata)./sqrt(numel(ydata));
        ny(xi) = length(ydata);
    end
%% create a structure for this curve, add errorbar plot to the current axes
    curvestructarray(zi) = ...
        struct('x', xs{zi}, 'ymean', ymean, 'ysem', ysem, 'nmeasures', ny, 'z', zs(zi), 'plotsym', plotsym{mod(zi-1, length(plotsym))+1}  );
    errorbar(curvestructarray(zi).x, curvestructarray(zi).ymean, curvestructarray(zi).ysem,...
        curvestructarray(zi).plotsym ,'MarkerEdgeColor',curvestructarray(zi).plotsym(2), 'MarkerFaceColor', curvestructarray(zi).plotsym (2),'MarkerSize',10);
    hold on
end
if ~isempty(paramcolindex)  % want a legend (a Z column was specified)
    legend(legs);
end
if exist('datasourcename','var')    % Case of file input
    title(['reduced data from file ', datasourcename]);
end
xlabel(['column ',num2str(xcolindex),' value']);
ylabel(['column ',num2str(ycolindex),' value']);
shg % display graph
grid