function [FMat,TMat,PFNMat,PTNMat,PFPMat,PTPMat,RMat,BatchFilePath,NChanExtra]=RepCon(InFileMat,NFilesPerCond,Weights,PermStatus,...
    MinPoint,MaxPoint,CalcBaseStatus,MinBase,MaxBase,TempFiltStruc,SpatFiltStruc,CalcGradStatus)
%   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.                                                                      
%==================================================================
%
%	RepCon.m
%
%	Nathan Weisz & Markus Junghoefer [2004]
%
%
%	Function definition
%	
%
%
%
%=================================================================
if nargin<12; CalcGradStatus=[];end
if nargin<11; SpatFiltStruc=[];end
if nargin<10; TempFiltStruc=[];end
if nargin<9; MaxBase=[]; end
if nargin<8; MinBase=[]; end
if nargin<7; CalcBaseStatus=[]; end
if nargin<6; MaxPoint=[]; end
if nargin<5; MinPoint=[]; end
if nargin<4; PermStatus=[]; end
if nargin<3; Weights=[]; end
if nargin<2; NFilesPerCond=[]; end
if nargin<1; InFileMat=[]; end


if isempty(TempFiltStruc); TempFiltStatus=0; else; TempFiltStatus=1; end
if isempty(SpatFiltStruc); SpatFiltStatus=0; else; SpatFiltStatus=1; end
if isempty(CalcGradStatus); CalcGradStatus=0; end

[PermStatus]=IfEmptyInputBo('Do you want to calculate','the permutation F-distribution ?',PermStatus,1,0);

[CalcBaseStatus]=IfEmptyInputBo('Do you want to subtract a baseline ?',[],CalcBaseStatus,1,1);
if CalcBaseStatus & (isempty(MinBase) | isempty(MaxBase))
    fprintf('\n')
    [NPointsBaseStatus]=IfEmptyInputBo('Do you want to use','the whole interval as baseline ?',[],1,0);
    if NPointsBaseStatus; 
        MinBase=[]; MaxBase=[]; 
    else
        [MinBase]=IfEmptyInputVal('Please insert the baseline start point:',[],MinBase,1,0);
        [MaxBase]=IfEmptyInputVal('Please insert the baseline start point:',[],MaxBase,1,0);
    end
end

[DefFilePath]=SetDefPath('r','*.rep*');
[DefFilePath]=SwitchGrandBatch(DefFilePath,'Grand2Batch',1);
StrTmp='Choose avg files or batch file:';

[NFilesIn,InFileMat,NonUse,BatchFilePath]=ReadFileNames(InFileMat,DefFilePath,StrTmp);
if NFilesIn==0; FMat=[];TMat=[];PFNMat=[];PTNMat=[];PFPMat=[];PTPMat=[];RMat=[];BatchFilePath=[];NChanExtra=[];return; end
[BatchFile,BatchPath]=SepFilePath(BatchFilePath);
SetDefPath(2,BatchPath);

NCond=length(Weights);
if isempty(NFilesPerCond)
    NFilesPerCond=NFilesIn./NCond;
    DefStr=num2str(NFilesPerCond);
    for i=2:NCond
        DefStr=[DefStr,' ',num2str(NFilesPerCond)];
    end
    Str=inputdlg({'Number of files per condition:'},'Insert vector:',1,{DefStr});
    NFilesPerCond=round(str2num(char(Str)));
end

Message=char('');
Message=char(Message,'Read data from files ...');
Message=char(Message,'');
Message=char(Message,'');
if CalcBaseStatus
    Message=char(Message,'Calculate baseline ...');
else
    Message=char(Message,'No baseline subtraction ...');
end
if TempFiltStatus 
    Message=char(Message,'Temporal filtering ...');
else
    Message=char(Message,'No temporal filtering ...');
end
if SpatFiltStatus
    Message=char(Message,'Spatial filtering ...');
else
    Message=char(Message,'No spatial filtering ...');
end
if CalcGradStatus
    Message=char(Message,'Calculate gradient ...');
else
    Message=char(Message,'No gradient ...');
end
hMsgbox1=msgbox(Message,'Info:','help');
hChildMsgbox1=get(hMsgbox1,'Children');
set(hChildMsgbox1(end),'Visible','off');
if ~strcmp(computer,'MAC2')
    hChild2Msgbox1=get(hChildMsgbox1(2),'Children');
end
pause(.1);

Message=[];
ActFileIndex=0;
for CondInd=1:NCond
    CondStr=int2str(CondInd);
    WeightsStr=int2str(Weights(CondInd));
    if CondInd<10; CondStr=['0',CondStr]; end
    Message=char(Message,['Condition: ',CondStr,'; Weight: ',WeightsStr]);
    Message=char(Message,'');
    for FileInd=1:NFilesPerCond(CondInd)
        ActFileIndex=ActFileIndex+1;
        [InFile,InPath,InFilePath]=GetFileNameOfMat(InFileMat,ActFileIndex);
        FileIndStr=int2str(FileInd);
        if FileInd<10; FileIndStr=['0',FileIndStr]; end
        Message=char(Message,[FileIndStr,': ',InFile]);
        [AvgMat,File,Path,FilePath,NTrialAvgVec,StdMat,...
                SampRate,AvgRef,Version,MedMedRawVec,MedMedAvgVec,...
                EegMegStatus,NChanExtra]=ReadAvgFile(InFilePath);
        [NChan,NPoints]=size(AvgMat);
        if NChanExtra>0
            AvgMat=AvgMat(1:NChan-NChanExtra,:);
            NChan=NChan-NChanExtra;
        end
        if CalcBaseStatus
            if isempty(MinBase); MinBase=1; end
            if isempty(MaxBase); MaxBase=NPoints1; end
            fprintf('Subtract baseline interval of points %g - %g\n',MinBase,MaxBase);
            [AvgMat]=CalcBaseline(AvgMat,MinBase,MaxBase,CalcBaseStatus);
        else
            fprintf('No baseline subtraction ...\n');
        end
        if TempFiltStatus 
            LowFiltCoeff=TempFiltStruc.LowFiltCoeff;
            HighFiltCoeff=TempFiltStruc.HighFiltCoeff;
            if ~isempty(LowFiltCoeff)
                fprintf('Temporal lowpass filtering ...\n');
                for i=1:NChan
                    AvgMat(i,:)=(filtfilt(LowFiltCoeff(1,:),LowFiltCoeff(2,:),AvgMat(i,:)'))';
                end
            end
            if ~isempty(HighFiltCoeff)
                fprintf('Temporal highpass filtering ...\n');
                for i=1:NChan
                    AvgMat(i,:)=(filtfilt(HighFiltCoeff(1,:),HighFiltCoeff(2,:),AvgMat(i,:)'))';
                end
            end
        else
            fprintf('No temporal filtering ...\n');
        end
        if SpatFiltStatus
            fprintf('Apply spatial lowpass filter. Sigma: %g\n',SpatFiltStruc.Sigma./pi.*180);
            AvgMat=SpatFilt(AvgMat,SpatFiltStruc.PhiMat,SpatFiltStruc.Sigma,1);
        else
            fprintf('No spatial filtering ...\n');
        end
        if CalcGradStatus
            fprintf('Calculate gradient\n');
            AvgMat=gradient(AvgMat);
        else
            fprintf('No gradient ...\n');
        end
        if CondInd==1 & FileInd==1; 
            AllAvgMat=zeros(max(NFilesPerCond),NCond,NChan,NPoints);
        end
        AllAvgMat(FileInd,CondInd,:,:)=AvgMat;
    end
    if CondInd<NCond; Message=char(Message,''); end
end
MessageSize=size(Message,1);
if MessageSize>160
    NewMessage=[];
    if MessageSize./3~=floor(MessageSize./3)
        Message=char(Message,'');
        MessageSize=MessageSize+1;
        if MessageSize./3~=floor(MessageSize./3)
            Message=char(Message,'');
            MessageSize=MessageSize+1;
        end
    end
    NewMessageSize=MessageSize./3;
    for i=1:NewMessageSize
        NewMessage=char(NewMessage,[Message(i,:),'       ;       ',Message(NewMessageSize+i,:),'       ;       ',Message(2.*NewMessageSize+i,:)]);
    end
    Message=NewMessage;
elseif MessageSize>80 & MessageSize<=160
    NewMessage=[];
    if MessageSize./2~=floor(MessageSize./2)
        Message=char(Message,'');
        MessageSize=MessageSize+1;
    end
    NewMessageSize=MessageSize./2;
    for i=1:NewMessageSize
        NewMessage=char(NewMessage,[Message(i,:),'       ;       ',Message(NewMessageSize+i,:)]);
    end
    Message=NewMessage;
end

hMsgbox=msgbox(Message,'Info:','help');
set(hMsgbox,'Units','Normal');
MsgBoxPos=get(hMsgbox,'Position');
MsgBoxPos(1:2)=[0.1 0];
set(hMsgbox,'Position',MsgBoxPos);
pause(.1)
FMat=zeros(NChan,NPoints);
TMat=zeros(NChan,NPoints);
PFNMat=zeros(NChan,NPoints);
PTNMat=zeros(NChan,NPoints);
PFPMat=zeros(NChan,NPoints);
PTPMat=zeros(NChan,NPoints);
RMat=zeros(NChan,NPoints);

Message=char('');
Message=char(Message,'Points done ...');
Message=char(Message,'');
Message=char(Message,'');

if isempty(MinPoint); MinPoint=1; end
if isempty(MaxPoint); MaxPoint=NPoints; end
for PointInd=MinPoint:MaxPoint
    if floor(PointInd./10)==PointInd./10
        fprintf(1,'%g of %g points done ...\n',PointInd-MinPoint+1,MaxPoint-MinPoint+1)
        MessageTmp=char(Message,[int2str(PointInd-MinPoint+1),' of ',int2str(MaxPoint-MinPoint+1)]);
        if strcmp(computer,'MAC2')
            set(hChildMsgbox1(2),'String',MessageTmp);
        else
            set(hChild2Msgbox1,'String',MessageTmp);
        end
        pause(.1);
    end
    RepMat=squeeze(AllAvgMat(:,:,:,PointInd));
    [NSub,NCond,NChan]=size(RepMat);
    
    for ChanInd=1:NChan     
        %Berechnung von Mean square contrast
        SumRepMat = sum(RepMat(:,:,ChanInd));%summierte ueber Vps 
        SumWeights2=sum(Weights.^2);
        L=sum(SumRepMat.*Weights);
        L2=L.^2;     
        MScont=L2/(NSub.*SumWeights2);
        
        %Corrections of data matrix   
        %for Grand Mean     
        corGM=RepMat(:,:,ChanInd)-mean(mean(RepMat(:,:,ChanInd)));
        
        %for row effect    
        rowef=mean(corGM,2);   
        [NSub,NCond]=size(RepMat(:,:,ChanInd));  
        for CondInd=1:NCond
            corGMR(1:NSub,CondInd)=corGM(1:NSub,CondInd)-rowef;
        end  
        
        %for column effect 
        colef = mean(corGMR);
        for SubInd =1:NSub
            corGMRC(SubInd,1:NCond)=corGMR(SubInd,1:NCond)-colef;
        end
        
        %sum of squares for condition x subjects
        %this is used as error term for calculation of F. one also could
        %calculate seperate error terms for each trend type. see Rosnow &
        %Rosenthal for details. the approach here is sufficient in most cases.
        
        SScs=sum(sum(corGMRC.^2));
        dfcs=NSub.*(NCond-1);
        MScs=SScs./dfcs;
        Fcont=MScont./MScs;
        Tcont=L/sqrt(SScs.*SumWeights2);
        FMat(ChanInd,PointInd)=Fcont; 
        PFNMat(ChanInd,PointInd)=fcdf(Fcont,1,dfcs);
        TMat(ChanInd,PointInd)=Tcont; 
        if Tcont>=0
            Pcont=tcdf(Tcont,dfcs);
        else
            Pcont=-tcdf(-Tcont,dfcs);
        end
        PTNMat(ChanInd,PointInd)=Pcont;
        %Calculat effect size r(contrast). Use of df of pooled error term will lead
        %to underestimation of effect size. df has to be divided by factor number
        %of conditions minus one.
        RMat(ChanInd,PointInd)=sqrt(Fcont/(Fcont+(dfcs/(NCond-1))));      
    end
end
MeanF=0;
pF=0;
while pF<.95
    MeanF=MeanF+.01;
    pF=fcdf(MeanF,1,dfcs);
end
MeanT=0;
pT=0;
while pT<.95
    MeanT=MeanT+.01;
    pT=tcdf(MeanT,dfcs);
end
MessageTmp=char(Message,'')
MessageTmp=char(MessageTmp,['95% F value: ',num2str(MeanF)]);
MessageTmp=char(MessageTmp,['95% t value: ',num2str(MeanT)]);
if strcmp(computer,'MAC2')
    set(hChildMsgbox1(2),'String',MessageTmp);
else
    set(hChild2Msgbox1,'String',MessageTmp);
end
pause(.1);

if PermStatus
    PermVec=randperm(NChan.*NPoints);
    MinPerm=100;
    MaxPerm=200;
    MaxCondPerm=1000;
    FPermMat=zeros(MaxPerm,MaxCondPerm);
    FPermVec=zeros(MaxPerm,1);
    MeanFPermVec=zeros(MaxPerm,1);
    DiffMeanFPerm=1;
    PermInd=0;
    
    hPermFig=figure('Units','normal','Position',[.66 .66 .33 .33],'Name','Permutation','Color',[.99 .99 .99],'NumberTitle','off');
    while PermInd<=MinPerm | (DiffMeanFPerm>0.001 & PermInd<MaxPerm) 
        PermInd=PermInd+1;
        ChanInd=ceil(PermVec(PermInd)./NPoints)
        PointInd=PermVec(PermInd)-(ChanInd-1).*NPoints
        
        RepMat=squeeze(AllAvgMat(:,:,ChanInd,PointInd));
        RandRepMat=RepMat;
        
        for CondPermInd=1:MaxCondPerm
            for SubInd=1:NSub
                RandRepMat(SubInd,:)=RepMat(SubInd,randperm(NCond));
            end
            RepMat=RandRepMat;
            %Berechnung von Mean square contrast
            SumRepMat = sum(RepMat);%summiert ueber Vps 
            SumWeights2=sum(Weights.^2);
            L=sum(SumRepMat.*Weights);
            L2=L.^2;     
            MScont = L2/(NSub*SumWeights2);
            
            %Corrections of data matrix for Grand Mean     
            corGM=RepMat-mean(mean(RepMat));
            
            %for row effect    
            rowef=mean(corGM,2);   
            for CondInd=1:NCond
                corGMR(1:NSub,CondInd) = corGM(1:NSub,CondInd)-rowef;
            end  
            %for column effect 
            colef = mean(corGMR);
            for SubInd =1:NSub
                corGMRC(SubInd,1:NCond)=corGMR(SubInd,1:NCond)-colef;
            end
            
            %sum of squares for condition x subjects
            %this is used as error term for calculation of F. one also could
            %calculate seperate error terms for each trend type. see Rosnow &
            %Rosenthal for details. the approach here is sufficient in most cases.
            
            SScs = sum(sum(corGMRC.^2));
            dfcs = NSub.*(NCond-1);
            MScs = SScs./dfcs;
            Fcont = MScont/MScs;
            TCont = L./sqrt(SScs.*SumWeights2);
            FCondPermVec(CondPermInd)=Fcont;   
            TCondPermVec(CondPermInd)=TCont; 
        end
        SortFCondPermVec=sort(FCondPermVec);
        SortTCondPermVec=sort(TCondPermVec);
        FPermMat(PermInd,:)=SortFCondPermVec;
        TPermMat(PermInd,:)=SortTCondPermVec;
        FPermVec(PermInd)=SortFCondPermVec(floor(0.95.*MaxCondPerm));
        TPermVec(PermInd)=SortTCondPermVec(floor(0.95.*MaxCondPerm));
        MeanFPermVec(PermInd)=mean(FPermVec(1:PermInd));
        MeanAbsTPermVec(PermInd)=mean(abs(TPermVec(1:PermInd)));
        plot(MeanFPermVec(1:PermInd))
        hold on
        title('Mean 95% F-value')
        pause(.1)
        if PermInd>1
            DiffMeanFPerm=abs(MeanFPermVec(PermInd)-MeanFPermVec(PermInd-1));
            fprintf(1,'%g permutations done ...\n',PermInd)
            fprintf(1,'Mean F: %g\n',MeanFPermVec(PermInd))
            fprintf(1,'Mean F diff: %g\n',DiffMeanFPerm)
            
            MessageTmp=char(Message,   ['Permutations: ',int2str(PermInd)]);
            MessageTmp=char(MessageTmp,['Mean 95% F      : ',num2str(MeanFPermVec(PermInd))]);
            MessageTmp=char(MessageTmp,['Mean F diff  : ',num2str(DiffMeanFPerm)]);
            if strcmp(computer,'MAC2')
                set(hChildMsgbox(2),'String',MessageTmp);
            else
                set(hChild2Msgbox,'String',MessageTmp);
            end
            pause(.1);
        end       
    end
    FPermMat=FPermMat(1:PermInd,:);
    MeanFPermMat=mean(FPermMat);
    MeanTPermMat=mean(TPermMat);
    plot(MeanFPermMat)
    hold on
    title('F - Distribution')
    PFPMat=zeros(NChan,NPoints);
    for CondPermInd=2:MaxCondPerm-1 
        Tmp=find(FMat>=mean(MeanFPermMat(CondPermInd-1:CondPermInd)) & FMat<mean(MeanFPermMat(CondPermInd:CondPermInd+1))); 
        if ~isempty(Tmp)
            PFPMat(Tmp)=CondPermInd./MaxCondPerm.*ones(size(Tmp));
        end
    end
    Tmp=find(FMat>MeanFPermMat(MaxCondPerm));
    if ~isempty(Tmp)
        PFPMat(Tmp)=ones(size(Tmp));
    end
    CenterInd=MaxCondPerm./2; %500
    for CondPermInd=CenterInd+2:MaxCondPerm-1 %502:999
        Tmp=find(TMat>=mean(MeanTPermMat(CondPermInd-1:CondPermInd)) & TMat<mean(MeanTPermMat(CondPermInd:CondPermInd+1))); 
        if ~isempty(Tmp)
            PTPMat(Tmp)=(CondPermInd-CenterInd)./(MaxCondPerm-CenterInd).*ones(size(Tmp));
        end
    end
    Tmp=find(TMat>MeanTPermMat(MaxCondPerm));
    if ~isempty(Tmp)
        PTPMat(Tmp)=ones(size(Tmp));
    end
    for CondPermInd=CenterInd-1:-1:2 %499:2
        Tmp=find(TMat<=mean(MeanTPermMat(CondPermInd:CondPermInd+1)) & TMat>mean(MeanTPermMat(CondPermInd:CondPermInd-1))); 
        if ~isempty(Tmp)
            PTPMat(Tmp)=-1.*(CenterInd-CondPermInd)./(MaxCondPerm-CenterInd).*ones(size(Tmp));
        end
    end
    Tmp=find(TMat<MeanTPermMat(1));
    if ~isempty(Tmp)
        PTPMat(Tmp)=ones(size(Tmp));
    end
    MeanF=MeanFPermVec(PermInd);
    MeanT=MeanAbsTPermVec(PermInd);
else
    PFPMat=zeros(NChan,NPoints); 
    PTPMat=zeros(NChan,NPoints); 
end
MessageTmp=char(Message,'')
MessageTmp=char(Message,['95% F value: ',num2str(MeanF)]);
MessageTmp=char(Message,['95% t value: ',num2str(MeanT)]);
if strcmp(computer,'MAC2')
    set(hChildMsgbox1(2),'String',MessageTmp);
else
    set(hChild2Msgbox1,'String',MessageTmp);
end
pause(.1);
sprintf('concerning F: degrees of freedom are 1,%g',dfcs);
sprintf('concerning t: degrees of freedom is %g',dfcs);
sprintf('concerning p: 95% F is %g',MeanF);
sprintf('concerning p: 95% t is %g',MeanT);
return;
