function [set]=CRCM1D(dire,n1,n2,n3)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%R.Magalhaes, Paris, July 2019
%Function to build the cross-realization correlation matrix
%as described in Yang et al (2007)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%This function calculates the pearson correlation
%coefficient between all pairs of components for the n
%permuations of the maps present in the directory
%building the Cross Realization Correlation Matrix (CRCM)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Must choose the range of values to build correlation matrix
%and step size (e.g between 5 and 400 permutations in steps of 5
%However, likely can be run only once for the higher permutation value as
%all other steps will be subcases of that one
%Usage: CRCM('directory',min number of permutations, size step,max number of permutations);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%loaded file must have a structure called "data" with one entry per
%permutation with file and name subfields. Name contains the id of the
%permutation and file the 4D output from ICA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%This software is provided under a Free for Non-Commercial Use Only Agreement and a Non-Commercial use agreement
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%contact: ricardo.magalhaes@cea.fr


load(strcat(dire,'/CRCM_data401'),'data');
for n=n1:n2:n3
    index=[];
    counti=1;
    countj=1;
    for i=1:1
        i
        data1=data(i).file;
        for j=1:n
            j
            data2=data(j).file;
            c=calc_correlations(data1,data2);
            if(i==j); c=zeros(size(c));end
            set(i,j).mat=c;
            set(i,j).name=strcat(int2str(i),'_',int2str(j));
            set(j,i).mat=c';
            set(j,i).name=strcat(int2str(j),'_',int2str(i));
            
        end
    end
    for i=1:size(set,1)
        sizes(i)=size(set(1,i).mat,2);
        index=[index [1:sizes(i);ones(1,sizes(i))*i]];
    end
    crcm=zeros(sum(sizes),sum(sizes));
          counti=1;
          for i=1:size(set,1)
              countj=1;
              for j=1:size(set,2)
                  crcm(counti:(counti+sizes(i)-1),countj:(countj+sizes(j))-1)=set(i,j).mat;
                  countj=countj+sizes(j);
              end
              counti=counti+sizes(i);
          end
    for x=1:1
        for y=1:length(set)
            set(x,y).indexx=[1:size(set(x,y).mat,1)];
            set(x,y).indexy=[1:size(set(x,y).mat,2)];
        end
    end
    save(strcat('CRCM_mat1D',int2str(n)),'set','sizes','index','crcm','-v7.3');
end
end

function mat=calc_correlations(set1,set2)
for i=1:size(set1,4)
    for j=1:size(set2,4)
        c=corrcoef(squeeze(set1(:,:,:,i)),squeeze(set2(:,:,:,j)));
        mat(i,j)=c(1,2);
    end
end
end