function [interactiondata]=calc_nway_interaction(anovadata,factor,n)

%    calc_nway_interaction
%
%   'calc_nway_interaction(anovadata,nrofsubjects,validfactor,n)'
%   calculates n-way-interaction for the data given in anovadata.
%   the number of rows in anovadata corresponds the number of 
%   subjects. 'factor' is a cell containing the name of the factors
%   in the first column and the number of gradations in the second
%   column. 'n' specifies the order of the interactions to be 
%   calculated. Factors must be in hierarchical order, starting 
%   from the highest level.  


%   This software is protected by german copyright and international treaties.             
%   Copyright 2004 Markus Junghöfer & Peter Peyk. All Rights Reserved.                     
%                                                                                          
%   THIS SOFTWARE AND DOCUMENTATION IS PROVIDED "AS IS," AND COPYRIGHT HOLDERS MAKE        
%   NO REPRESENTATIONS OR WARRANTIES, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO,    
%   WARRANTIES OF MERCHANTABILITY OR FITNESS FOR ANY PARTICULAR PURPOSE OR THAT THE USE    
%   OF THE SOFTWARE OR DOCUMENTATION WILL NOT INFRINGE ANY THIRD PARTY PATENTS, COPYRIGHTS,
%   TRADEMARKS OR OTHER RIGHTS. COPYRIGHT HOLDERS WILL NOT BE LIABLE FOR ANY DIRECT,       
%   INDIRECT, SPECIAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF ANY USE OF THE SOFTWARE      
%   OR DOCUMENTATION.                                                                      



nrofsubjects = size(anovadata,1);
nroffactors = size(factor,1);
nrofinteractions = nchoosek(size(factor,1),n);
GrandMean = mean(mean(anovadata));
CellMean = mean(anovadata,1);

%analyze factor structure
%------------------------------
totcells = 1;
for j=1:size(factor,1)
    totcells = totcells * factor{j,2};
end
if size(anovadata,2)~= totcells
    uiwait(errordlg('Number of columns of data does not match the product of factor gradations!','Error calculating factor interaction!:'));
    interactiondata = [];
    return
end

factorstructure = zeros(nroffactors,totcells);
for factorzaehler=nroffactors:-1:1
    
    %check gradations
    nrofgradations = factor{factorzaehler,2};
    if round(totcells/nrofgradations)~=totcells/nrofgradations
        uiwait(errordlg(['Number of cells cannot be divided in n-gradations for factor ',factor{factorzaehler,1},'!'],'Error:'));
    end
    %step size
    step = 1;
    for bottomzaehler = factorzaehler+1:nroffactors
        step = step * factor{bottomzaehler,2};
    end
    %update factorstructure
    gradationsvec = linspace(1,nrofgradations,nrofgradations);
    gradationzaehler = 1;
    for columnzaehler=1:step:totcells-(step-1)
        factorstructure(factorzaehler,columnzaehler:columnzaehler+(step-1))=gradationsvec(gradationzaehler);
        gradationzaehler = gradationzaehler +1;
        if gradationzaehler > nrofgradations
            gradationzaehler = 1;
        end
    end
end
 
%types of interactions
combinations = nchoosek(linspace(1,nroffactors,nroffactors),n);

%----------------------------------------------------------------------------
%interactiondata:  1: condition   2:QS   3:df   4: errorQS  5:errordf  6:F-Value:   7:Significance;
%----------------------------------------------------------------------------
interactiondata = cell(size(combinations,1),6);

for combinationzaehler=1:size(combinations,1)
    interactiondata{combinationzaehler,1}=combinations(combinationzaehler,:);
    
    factorstructure_ia = factorstructure;
    komponenten = cell(n,2);
    nrofinteractioncells = 1;
    sizestring = '';
    conditionindexstring = '';
    for komponentenzaehler=1:n
        %Name der Komponente
        komponenten{komponentenzaehler,1}=factor{combinations(combinationzaehler,komponentenzaehler),1};
        %nrofgradations der Komponente
        komponenten{komponentenzaehler,2}=factor{combinations(combinationzaehler,komponentenzaehler),2};
        if ~isempty(sizestring)
            sizestring = [sizestring,',',num2str(komponenten{komponentenzaehler,2})];
            conditionindexstring = [conditionindexstring,',','factorstructure_ia(',num2str(komponentenzaehler),',r)'];
        else
            sizestring = [num2str(komponenten{komponentenzaehler,2})];
            conditionindexstring = ['factorstructure_ia(',num2str(komponentenzaehler),',r)'];
        end
        nrofinteractioncells = nrofinteractioncells * komponenten{komponentenzaehler,2};
    end
    
    evalstring = ['InteractionMean = zeros(',sizestring,');']; %--->InteractionMean = zeros(.....);
    eval(evalstring);
    
    %conditionmat
    %--------------
    conditionmat = zeros(size(InteractionMean));
    conditionmat(1:end)=linspace(1,nrofinteractioncells,nrofinteractioncells);
    
    %conditionvec
    %---------------
    conditionvec = zeros(1,size(anovadata,2));
    %restzeilen löschen
    deletevec=[];
    deletezaehler=1;
    for r=1:nroffactors
       if isempty(find(combinations(combinationzaehler,:)==r))
           deletevec(deletezaehler)=r;
           deletezaehler=deletezaehler+1;
       end
    end
    factorstructure_ia(deletevec,:)=[];
    evalstring = ['conditionvec(1,r)=conditionmat(',conditionindexstring,');'];
    for r=1:size(factorstructure_ia,2)
        eval(evalstring);
    end
    
    %calculate conditionmeans, QS & df |  calculate test-qs
    %-------------------------------------------------------
    ConditionQS = 0;
    TestQS = 0;
    for r=1:nrofinteractioncells
        InteractionMean(r)=mean(CellMean(find(conditionvec==r)));
        ConditionQS = ConditionQS + ((InteractionMean(r)-GrandMean)^2);
        for s=1:nrofsubjects
            TestQS = TestQS + (mean(anovadata(s,find(conditionvec==r)))-GrandMean)^2;
        end
    end
    interactiondata{combinationzaehler,2}=ConditionQS;
    interactiondata{combinationzaehler,3}=prod(size(InteractionMean)-1);%df
    conditiondf = interactiondata{combinationzaehler,3};
    interactiondata{combinationzaehler,4}=TestQS;%error
    interactiondata{combinationzaehler,5}=prod([(size(InteractionMean)-1),nrofsubjects-1]);%errordf
    testdf =  interactiondata{combinationzaehler,5}; 
    interactiondata{combinationzaehler,6}=(ConditionQS/conditiondf)/(TestQS/testdf);% F-Value
    f_value = interactiondata{combinationzaehler,6};
    interactiondata{combinationzaehler,7}=fcdf(f_value,conditiondf,testdf); %p
    
      
    
end
    
    
        



return