%  ANOVA2_CELL - compute F-values in cell array using ANOVA.
%
% Usage:
%    >> [FC FR FI dfc dfr dfi] = anova2_cell( data );
%
% Inputs:
%   data       = data consisting of PAIRED arrays to be compared. The last 
%                dimension of the data array is used to compute ANOVA.
% Outputs:
%   FC   - F-value for columns.
%   FR   - F-value for rows.
%   FI   - F-value for interaction.
%   dfc  - degree of freedom for columns.
%   dfr  - degree of freedom for rows.
%   dfi  - degree of freedom for interaction.
%
% Note: the advantage over the ANOVA2 function of Matlab statistical
%       toolbox is that this function works on arrays (see examples). Note
%       also that you still need the statistical toolbox to assess
%       significance using the FCDF function. The other advantage is that
%       this function will work with complex numbers.
%
% Example:
%   a = { rand(1,10) rand(1,10) rand(1,10); rand(1,10) rand(1,10) rand(1,10) }
%   [FC FR FI dfc dfr dfi] = anova2_cell(a)
%   signifC = 1-fcdf(FC, dfc(1), dfc(2))
%   signifR = 1-fcdf(FR, dfr(1), dfr(2))
%   signifI = 1-fcdf(FI, dfi(1), dfi(2))
%
%   % for comparison 
%   anova2(  [ a{1,1}' a{1,2}' a{1,3}'; a{2,1}' a{2,2}' a{2,3}' ], 10) 
%
%   b = { [ a{1,1}; a{1,1} ] [ a{1,2}; a{1,2} ] [ a{1,3}; a{1,3} ];
%         [ a{2,1}; a{2,1} ] [ a{2,2}; a{2,2} ] [ a{2,3}; a{2,3} ] }
%   [FC FR FI dfc dfr dfi] = anova2_cell(b)
%
%   c{1,1} = reshape(repmat(b{1,1}, [2 1]),2,2,10);
%   c{1,2} = reshape(repmat(b{1,2}, [2 1]),2,2,10);
%   c{1,3} = reshape(repmat(b{1,3}, [2 1]),2,2,10);
%   c{2,3} = reshape(repmat(b{2,3}, [2 1]),2,2,10);
%   c{2,2} = reshape(repmat(b{2,2}, [2 1]),2,2,10);
%   c{2,1} = reshape(repmat(b{2,1}, [2 1]),2,2,10)
%   [FC FR FI dfc dfr dfi] = anova2_cell(c)
%
% Author: Arnaud Delorme, SCCN/INC/UCSD, La Jolla, 2005
%
% Reference:
%   Schaum's outlines in statistics (3rd edition). 1999. Mc Graw-Hill.

% Copyright (C) Arnaud Delorme
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

function [FC, FR, FI, freeC, freeR, freeI] = anova2_cell(data)
    
    % compute all means and all std
    % -----------------------------
    a = size(data,1);
    b = size(data,2);
    nd = myndims( data{1} );
    c  = size(data{1}, nd);
    sz = size(data{1});
    
    % dataabs if for complex data only
    % --------------------------------
    dataabs = data;
    if ~isreal(data{1})
        for i = 1:a
            for ii = 1:b
                dataabs{i,ii} = abs(data{i,ii});
            end
        end
    end
    
    VE = zeros( [ sz(1:end-1) 1], 'single' );
    m  = zeros( [ sz(1:end-1) size(data) ], 'single' );
    for i = 1:a
        for ii = 1:b
            tmpm = mymean(data{i,ii}, nd);
            switch nd
                case 1, m(i,ii) = tmpm;
                case 2, m(:,i,ii) = tmpm;
                case 3, m(:,:,i,ii) = tmpm;
                case 4, m(:,:,:,i,ii) = tmpm;
                case 5, m(:,:,:,:,i,ii) = tmpm;
                case 6, m(:,:,:,:,:,i,ii) = tmpm;
                case 7, m(:,:,:,:,:,:,i,ii) = tmpm;
                otherwise error('Dimension not supported');
            end
            VE = VE+sum( bsxfun(@minus, dataabs{i,ii}, tmpm).^2, nd);
        end
    end
    X  = mean(mean(m,nd+1),nd);
    Xj = mean(m,nd+1);
    Xk = mean(m,nd);
    VR = b*c*sum( bsxfun(@minus, Xj, X).^2, nd);
    VC = a*c*sum( bsxfun(@minus, Xk, X).^2, nd+1 );
    VI = c*sum( sum( bsxfun(@plus, bsxfun(@minus, bsxfun(@minus, m, Xj), Xk), X).^2, nd+1 ), nd );
        
% before bsxfun
%     VR = b*c*sum( (Xj-repmat(X, [ones(1,nd-1) size(Xj,nd  )])).^2, nd   );
%     VC = a*c*sum( (Xk-repmat(X, [ones(1,nd  ) size(Xk,nd+1)])).^2, nd+1 );
%     
%     Xj = repmat(Xj, [ones(1,nd)   size(m,nd+1) ]);
%     Xk = repmat(Xk, [ones(1,nd-1) size(m,nd)  1]);
%     VI = c*sum( sum( ( m - Xj - Xk + repmat(X, [ones(1,nd-1) size(m,nd) size(m,nd+1)]) ).^2, nd+1 ), nd );
    
%     if nd == 1
%         
%         VE = 0;
%         m  = zeros( size(data), 'single' );
%         for i = 1:a
%             for ii = 1:b
%                 m(i,ii) = mymean(data{i,ii});
%                 VE      = VE+sum( (dataabs{i,ii}-m(i,ii)).^2 );
%             end
%         end
%         X  = mean(mean(m));
%         Xj = mean(m,2);
%         Xk = mean(m,1);
%         VR = b*c*sum( (Xj-X).^2 );
%         VC = a*c*sum( (Xk-X).^2 );
%         
%         Xj = repmat(Xj, [1 size(m,2) ]);
%         Xk = repmat(Xk, [size(m,1)  1]);
%         VI = c*sum( sum( ( m - Xj - Xk + X ).^2 ) );
%         
%     elseif nd == 2
% 
%         VE = zeros( size(data{1},1),1, 'single');
%         m  = zeros( [ size(data{1},1) size(data) ], 'single' );
%         for i = 1:a
%             for ii = 1:b
%                 tmpm = mymean(data{i,ii}, 2);
%                 m(:,i,ii) = tmpm;
%                 VE        = VE+sum( (dataabs{i,ii}-repmat(tmpm, [1 size(data{i,ii},2)])).^2, 2);
%             end
%         end
%         X  = mean(mean(m,3),2);
%         Xj = mean(m,3);
%         Xk = mean(m,2);
%         VR = b*c*sum( (Xj-repmat(X, [1 size(Xj,2)])).^2, 2 );
%         VC = a*c*sum( (Xk-repmat(X, [1 1 size(Xk,3)])).^2, 3 );
%         
%         Xj = repmat(Xj, [1 1 size(m,3) ]);
%         Xk = repmat(Xk, [1 size(m,2)  1]);
%         VI = c*sum( sum( ( m - Xj - Xk + repmat(X, [1 size(m,2) size(m,3)]) ).^2, 3), 2 );
%         
%     elseif nd == 3
%         
%         VE = zeros( size(data{1},1), size(data{1},2), 'single' );
%         m  = zeros( [ size(data{1},1) size(data{1},2) size(data) ], 'single' );
%         for i = 1:a
%             for ii = 1:b
%                 tmpm = mymean(data{i,ii}, 3);
%                 m(:,:,i,ii) = tmpm;
%                 VE          = VE+sum( (dataabs{i,ii}-repmat(tmpm, [1 1 size(data{i,ii},3)])).^2, 3);
%             end
%         end
%         X  = mean(mean(m,4),3);
%         Xj = mean(m,4);
%         Xk = mean(m,3);
%         VR = b*c*sum( (Xj-repmat(X, [1 1 size(Xj,3)  ])).^2, 3 );
%         VC = a*c*sum( (Xk-repmat(X, [1 1 1 size(Xk,4)])).^2, 4 );
%         
%         Xj = repmat(Xj, [1 1 1 size(m,4) ]);
%         Xk = repmat(Xk, [1 1 size(m,3)  1]);
%         VI = c*sum( sum( ( m - Xj - Xk + repmat(X, [1 1 size(m,3) size(m,4)]) ).^2, 4 ), 3 );
%                 
%     else % nd == 4
%         
%         VE = zeros( size(data{1},1), size(data{1},2), size(data{1},3), 'single' );
%         m  = zeros( [ size(data{1},1) size(data{1},2) size(data{1},3) size(data) ], 'single' );
%         for i = 1:a
%             for ii = 1:b
%                 tmpm = mymean(data{i,ii}, 4);
%                 m(:,:,:,i,ii) = tmpm;
%                 VE            = VE+sum( (dataabs{i,ii}-repmat(tmpm, [1 1 1 size(data{i,ii},4)])).^2, 4);
%             end
%         end
%         X  = mean(mean(m,5),4);
%         Xj = mean(m,5);
%         Xk = mean(m,4);
%         VR = b*c*sum( (Xj-repmat(X, [1 1 1 size(Xj,4)  ])).^2, 4 );
%         VC = a*c*sum( (Xk-repmat(X, [1 1 1 1 size(Xk,5)])).^2, 5 );
%         
%         Xj = repmat(Xj, [1 1 1 1 size(m,5) ]);
%         Xk = repmat(Xk, [1 1 1 size(m,4)  1]);
%         VI = c*sum( sum( ( m - Xj - Xk + repmat(X, [1 1 1 size(m,4) size(m,5)]) ).^2, 5 ), 4 );
%                 
%     end
    
    SR2 = VR/(a-1);
    SC2 = VC/(b-1);
    SI2 = VI/(a-1)/(b-1);
    SE2 = VE/(a*b*(c-1));
    
    FR = SR2./SE2; % rows
    FC = SC2./SE2; % columns
    FI = SI2./SE2; % interaction
    
    freeR = [ a-1 a*b*(c-1) ];
    freeC = [ b-1 a*b*(c-1) ];
    freeI = [ (a-1)*(b-1) a*b*(c-1) ];

function val = myndims(a)
    if ndims(a) > 2
        val = ndims(a);
    else
        if size(a,1) == 1,
            val = 2;
        elseif size(a,2) == 1,
            val = 1;
        else
            val = 2;
        end
    end; 
  
function res = mymean( data, varargin) % deal with complex numbers
    res = mean( data, varargin{:});
    if ~isreal(data)
        res = abs( res );
    end
