function [datamatrix,data_labels,headers]=read_experimental_data(fname,vdims,cdims,selectors)
% function [datamatrix,data_labels,headers]=read_experimental_data(fname,vdims,cdims,selectors)
% MATLAB function to read experimental data gathered for 
% Ambisonics Book.
%
% fname is the file name
% vdims specifies the first two dimensions of the response, e.g. 1D or 2D
% cdims specifies the the condition parameter dimensions thereafter
%   so that vdims+cdims is the number of columns in the experimental data file
% selectors is a cell array of the length cdims
%   each entry can either be empty (all conditions will be selected in 
%   string-sorting order
%   or it can contain the selected or ordered condition strings
%
% Franz Zotter, 2018

fid=fopen(fname);
headers=textscan(fid,[repmat('%s ',1,vdims+cdims) '\n'],1);
data = textscan(fid, [repmat('%f ',1,vdims) repmat('%s ',1,cdims)]);
fclose(fid);

idx=true(length(data{1}),1);
for k=1:cdims
    if ~isempty(selectors{k})
        idxk=false(length(data{1}),1);
        selk=selectors{k};
        for kk=1:length(selk)
            idxk=idxk|strcmp(data{:,k+vdims},selk{kk});
        end
        idx=idx&idxk;
    end
end
for k=1:cdims+vdims
    datak=data{k}(idx);
    data{k}=datak;
end

data_labels=cell(1,cdims);
sz=zeros(1,cdims);
for k=1:cdims
    if isempty(selectors{k})
        data_labels{k}=unique(data{:,k+vdims});
    else
        data_labels{k}=selectors{k};
    end
    sz(k)=length(data_labels{k});
end


combinations=zeros(prod(sz),length(sz));
combinations(:,1)=0:prod(sz)-1;
for d=1:length(sz)-1
    combinations(:,d+1)=floor(combinations(:,d)/sz(d));
    combinations(:,d)=mod(combinations(:,d),sz(d));
end
combinations=combinations+1;
multiplicity=0;
for k=1:size(combinations,1)
    idx=true(length(data{1}),1);
    for d=1:cdims
        datalabelsk=data_labels{d}(combinations(k,d));
        idx=idx&strcmp(data{d+vdims},datalabelsk);
    end
    len=length(find(idx));
    multiplicity=max(len,multiplicity);
end


dims=[multiplicity vdims sz];
datamatrix=NaN(dims);


for k=1:size(combinations,1)
    idx=true(length(data{1}),1);
    for d=1:cdims
        datalabelsk=data_labels{d}(combinations(k,d));
        idx=idx&strcmp(data{d+vdims},datalabelsk);
    end    
    for dd=1:vdims
        v=data{dd}(idx);
        datamatrix(1:length(v),dd,k)=v;
    end
end