function [counts,vals] = xtab(varargin)

% XTAB -- Build (and print) a contingency table from data.
%
% [COUNTS,VALS] = XTAB(SET1,...,SETk)
%
% SET1,...,SETk are data matrices with the same number of rows (cases).
% Each set is treated as a separate dimension of the cross-tabulation.
% Each set can comprise of one or more columns (variables). When there
% are multiple variables in a set, the unique combinations thereof define
% the bins along this particular dimension.
%
% The output is a K-dimensional matrix COUNTS and, optionally, a cell
% array of K matrices of unique combinations for the corresponding data
% set.  COUNTS(I1,I2,...,Ik) contains the number of observations with
% with (composite) values [VALS{1}(I1),VALS{2}(I2),...,VALS{k}(Ik)].
%
% XTAB(...) without output arguments prints a cross-tabulation table.
% This option is available for up to two sets only (see XTAB1 and XTAB2).
%
% See also PRINT_XTAB, XTAB1, XTAB2, CHI2, CROSSTAB, UNIQUE.

% Original coding by Alex Petrov, Ohio State University
% $Revision: 1.1 $  $Date: 2002/01/18 17:40 $
%
% Part of the utils toolbox version 1.1 for MATLAB version 5 and up.
% http://alexpetrov.com/softw/utils/
% Copyright (c) Alexander Petrov 1999-2006, http://alexpetrov.com
% Please read the LICENSE and NO WARRANTY statement in ../utils_license.m

N_sets = length(varargin) ;

if (nargout==0)         % print requested ?
  if (N_sets==1)        % uni-dimensional data
    xtab1(varargin{1}) ;
    return
  elseif (N_sets==2)    % two-dimensional data
    xtab2(varargin{1},varargin{2}) ;
    return
  else                  % general case
    disp('PRINT_XTAB cannot handle more than two dimensions.');
  end
else  % export values and don't print anything
  if (N_sets==1)        % uni-dimensional data
    [counts,vals] = xtab1(varargin{1}) ;
    return
  elseif (N_sets==2)    % two-dimensional data
    [counts,vals] = xtab2(varargin{1},varargin{2}) ;
    return
  %else the rest of the file
  end
end

%%%%  General case: N_SETS >= 3  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%

N_cases = size(varargin{1},1) ;

%--  Determine unique values for each variable separately
vals = cell(1,N_sets) ;
keys = zeros(N_cases,N_sets) ;
N_vars = zeros(1,N_sets) ;    % number of columns within each data set
N_vals = zeros(1,N_sets) ;    % number of unique (composite) values
for k = 1:N_sets
   [N,N_vars(k)] = size(varargin{k}) ;
   if (N~=N_cases)
      error('There must be the same number of rows in each data set.') ;
   end
   if (N_vars(k)==1)   % simple values
      [vals{k},foo,keys(:,k)] = unique(varargin{k}) ;
   else        % multiple variables --> composite values
      [vals{k},foo,keys(:,k)] = unique(varargin{k},'rows') ;
   end
   N_vals(k) = size(vals{k},1) ;
end

%-- Intersect the (recoded) sets and count the unique tuples
[tuples,foo,t_keys] = unique(keys,'rows') ;
linear_counts = hist(t_keys,size(tuples,1))' ;

%-- Arrange LINEAR_COUNTS into a rectangular table
counts = zeros(N_vals) ;
tuples_aux = num2cell(tuples,1) ;       % convert to a cell array
idx = sub2ind(N_vals,tuples_aux{:}) ;   % linear index
counts(idx) = linear_counts ;

%% Return COUNTS and VALS
%%%%% End of XTAB.M
