function [mask] = permTest(EEG_dat,thresh,numPerm,MA,alpha)
% modified from Matt Nassar's code
if nargin<5
    alpha=.05;
end
load connectionMat
load chanlocs

plotfigure=0;  
          
% get clusters, cluster sizes, cluster masses for positive clusters (ie
% p<.05 in a one tailed positive Ctest):
posClusterInfo=getEEG_clusterSize_clean(EEG_dat, connectionMat, thresh, 'right');
% get clusters, cluster sizes, cluster masses for negative clusters (ie
% p<.05 in a one tailed negative test):
negClusterInfo=getEEG_clusterSize_clean(EEG_dat, connectionMat, thresh, 'left');

ts = posClusterInfo.tMap;

maxSize = [];
maxWt = [];
for i=1:numPerm
    % create a subject length array containing randomly assigned -1's and
    % 1s:
    permArray=ones(size(EEG_dat, 1), 1);
    permArray(logical(binornd(1, .5, size(EEG_dat, 1), 1)))=-1;
    
    % multiply each subject timeseries by the -1 or 1 assigned randomly on
    % this trial
    sz=size(EEG_dat);
    permMat=repmat(permArray, [1 sz(2:end)]);
    
    % get cluster statistics for permuted dataset:
    permClusterInfo=getEEG_clusterSize_clean(EEG_dat.*permMat, connectionMat, thresh, 'right');
    
    % store the maximum of the statistics, for use in null distribution:
    maxSize(i)=(max(permClusterInfo.clustSizeMap(:)));
    maxWt(i)=(max(permClusterInfo.clustWtMap(:)));
end

% For a two tailed test, find the the minimum cluster statistics necessary
% to beat 97.5% of the null distribution
sizeTh = prctile(maxSize,100*(1-alpha/2)) 
massTh = prctile(maxWt,100*(1-alpha/2)) 
if plotfigure
figure
subplot(1,2,1)
plot(sort(maxSize),'ok')
hold on
plot(.975*numPerm,sizeTh,'*r')
plot([1 numPerm],max(negClusterInfo.clustSizeMap(:))*[1 1],'b')
plot([1 numPerm],max(posClusterInfo.clustSizeMap(:)),'g')
hold off
subplot(1,2,2)
plot(sort(maxWt),'o')
hold on
plot(.975*numPerm,massTh,'*r')
plot([1 numPerm],max(negClusterInfo.clustWtMap(:))*[1 1],'b')
plot([1 numPerm],max(posClusterInfo.clustWtMap(:))*[1 1],'g')
hold off
end

% Look for clusters that beat the null:

% Based on "mass"  
% I like this statistic better!!!

gPos=unique(posClusterInfo.ID_map(posClusterInfo.clustWtMap	>massTh));
gNeg=unique(negClusterInfo.ID_map(negClusterInfo.clustWtMap	>massTh));

threshMap=zeros(size(negClusterInfo.tMap));
threshMap(posClusterInfo.clustWtMap	>massTh)=posClusterInfo.tMap(posClusterInfo.clustWtMap	>massTh);
threshMap(negClusterInfo.clustWtMap	>massTh)=negClusterInfo.tMap(negClusterInfo.clustWtMap	>massTh);

if plotfigure
% look at thresholded maps across time:
figure
imagesc(threshMap, [-5 5])

figure % look at all positive clusters
imagesc(posClusterInfo.ID_map)

figure % and just suprathreshold clusters
sigClust=posClusterInfo.ID_map;
sigClust(posClusterInfo.clustWtMap	<massTh)=-1;
imagesc(sigClust)
% close all
% 
figure % and just suprathreshold clusters
sigClust=negClusterInfo.ID_map;
sigClust(negClusterInfo.clustWtMap	<massTh)=-1;
imagesc(sigClust)
end





% We can also test significance based on "size". 
% I'm less into this stat, as it doesn't give you any credit for being
% "really signifant"... but we can look at it anyway:

sgPos=unique(posClusterInfo.ID_map(posClusterInfo.clustSizeMap	>sizeTh));
sgNeg=unique(negClusterInfo.ID_map(negClusterInfo.clustSizeMap	>sizeTh));

sthreshMap=zeros(size(negClusterInfo.tMap));
sthreshMap(posClusterInfo.clustSizeMap	>sizeTh)=posClusterInfo.tMap(posClusterInfo.clustSizeMap	>sizeTh);
sthreshMap(negClusterInfo.clustSizeMap	>sizeTh)=negClusterInfo.tMap(negClusterInfo.clustSizeMap	>sizeTh);




% Its sort of nice to look at the clusters in topoplots... which requires
% chanLocs... i was loading a separate EEGlab dataset just to get this, but
% i'm sure you guys have something like this lying around.
% load EEG_demoFile.mat
if plotfigure
a=figure;
timtopo(threshMap, chanlocs, 'title', ['mass thresh th =',num2str(thresh)])
b=figure
timtopo(sthreshMap, chanlocs, 'title', 'size thresh')
c=figure
timtopo(negClusterInfo.tMap, chanlocs, 'title', 'no thresh')
end
% you should be able to see the clusters in the topo maps... for anything that i've
% run the clusters look exactly the same using a mass threshold as they do
% using the size threshold.
t = 0;
mask = [];
allmask = 0*posClusterInfo.ID_map;
if plotfigure
figure
end
for i = gPos'
    t = t+1;
    mask{t} = (posClusterInfo.ID_map==i).*ts;
    if plotfigure
    subplot(length(gPos)+length(gNeg),1,t)
    imagesc(mask{t})
    end
    allmask = allmask+mask{t};
end
for i = gNeg'
    t = t+1;
    mask{t} = (negClusterInfo.ID_map==i).*ts;
    if plotfigure
    subplot(length(gPos)+length(gNeg),1,t)
    imagesc(mask{t})
    end
    allmask = allmask+mask{t};
end
mask{t+1} = allmask;
end