function [va vvr]=RAICAR1D_var(dire,n1,n2,n3)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%R.Magalhaes, Paris, July 2019
%Function to calculate the component order variance across the number of
%permutations (va) used as well as return a 4D matrix with the averaged matched
%components (vvr)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Must choose the range of values to build correlation matrix
%and step size (e.g between 5 and 400 permutations in steps of 5
%Usage: RAICARD1D_fin('directory',min number of permutations, size step,max number of permutations);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Loaded file can be created with matlab function RAICAR1D made available
%within the same set of functions 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%This software is provided under a Free for Non-Commercial Use Only Agreement and a Non-Commercial use agreement
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%contact: ricardo.magalhaes@cea.fr

c=1;
load(strcat(dire,'/CRCM_mat1D400'),'set','sizes','index','crcm');
load(strcat(dire,'/CRCM_data400'),'data');
setall=set;sizesall=sizes;indexall=index;crcmall=crcm;
clear set, sizes, index, crcm;
for i=n1:n2:n3
    i
    set=setall(1:i,1:i);
    sizes=sizesall(1:i);
    index=indexall(:,1:i);
    crcm=crcmall(1:sum(sizes),1:sum(sizes));
    load(strcat(dire,'/RAICAR1Dresultsnew',int2str(i)),'comp');
    vr=[];
    d=[];
    for ij=1:i
        for j=1:size(set(1,ij).mat,1)
            d=[d set(1,ij).mat(j,:)];
        end
    end
    hd=hist(abs(d),200);
    hdlog=log(hd);
    hdsm=smooth(hdlog,125,'sgolay');
    [m, in]=min(hdsm);
   
    for j=1:length(comp)
        j
        ck=1;
        if(comp(1,j).r(1,3)>in/200)
            vv=data(1,1).file(:,:,:,comp(1,j).r(1,1));
            vr(1,:,:,:)=vv;
            v(1,:)=vv(:)';
            ck=ck+1;
        end
        if(comp(1,j).r(2,3)>in/200)
            vv=data(1,comp(1,j).r(1,4)).file(:,:,:,comp(1,j).r(1,2));
            vr(2,:,:,:)=vv;
            v(2,:)=vv(:)';
            ck=ck+1;

        end
        
        for k=2:size(comp(1,j).r,1)
            if(comp(1,j).r(k,3)>in/200)
                vv=data(1,comp(1,j).r(k,4)).file(:,:,:,comp(1,j).r(k,2));
                v(ck,:)=vv(:)';
                vr(ck,:,:,:)=vv;
                ck=ck+1;
            end
        end
        if(isempty(vr));vr=zeros(1,48,48,39);end
        vvr(j,:,:,:)=mean(vr,1);
        vr=[];
        va1=var(v,1);
        va1(isnan(va1))=0;
        v=[];
        va2(j)=mean(va1);
        
    end
    vv2(c,:)=va2;
    va(c)=mean(va2);
    c=c+1;
end

end