% This script was developped by Jerome G. Prunier.
% Version 3.0 (the function can now be easily run in a loop)
  % This script computes the mean number of interpopulation individual exchanges per generation in CDpop simulations and provides a summary matrix in .txt format.
  % The script is based on the handling of individuals' ID in CDpop output grids.
  % Input parameters :
     % path : this is the directory you chose as an output folder when running CDpop and where the output matrix will be saved ; it should contain at least one subfolder named 'batchrun0mcrun0'),
     % batch : batchrun number
     % mc : mcrun number
  % Output matrix are named according to batchrun and mcrun numbers (e.g. 'batchrun0mcrun0_output.txt').

function CDpopinterpopExchange(path,batch,mc)
    folder1=path;
    batchrunnumber=batch;
    mcrunnumber=mc;
    outputname=['batchrun',num2str(batch),'mcrun',num2str(mc),'_output.txt'];
    folder2=[folder1,'\batchrun',num2str(batchrunnumber),'mcrun',num2str(mcrunnumber)];
    cd(folder2);
    path=[folder2,'\'];
    files = dir(folder2);
    fileIndex = find([files.isdir]==0);
    par=0;
    for i = 1:length(fileIndex)
        filename = files(fileIndex(i)).name;
        if sum(filename(1,1:4)=='grid')~=4
            continue
        else
            data = dataset('xlsfile',sprintf('%s%s', path,filename));
            tt=double(data.Subpopulation);
            prov=data.ID;
            if isa(prov,'double')==1
                continue
            else par=par+1;
                Key   = 'grid';
                Index = strfind(filename, Key);
                Generation = sscanf(filename(Index(1) + length(Key):end), '%g', 1);
                genlist(par)=Generation;
                disp(['          Compiling immigration at generation ',num2str(Generation)]);
                npop=size(unique(tt),1);
                for j=1:size(prov,1)
                    provstring=char(prov(j,1));
                    Key2   = 'op';
                    Index2 = strfind(provstring, Key2);
                    if sum(size(Index2))==0
                        continue
                    else
                        provpop = sscanf(provstring(Index2(1) + length(Key2):end), '%g', 1);
                        provvector(j,1,par)=tt(j);
                        provvector(j,2,par)=provpop;
                    end
                end
            end
        end
    end
    genlist=sort(genlist);
    provmatrix=zeros(npop,npop);
    for j=1:par
        if j==1
            temp=provvector(:,:,j);
        else temp2=provvector(:,:,j);
            temp=[temp;temp2];
        end
    end
    provtot=sortrows(temp,1);
    provtot(isnan(provtot(:,2))==1,:)=[];
    for j=1:npop
        provtottemp=provtot(provtot(:,1)==j,:);
        for k=1:size(provtottemp,1)
            provmatrix(j,provtottemp(k,2))=provmatrix(j,provtottemp(k,2))+1;
        end
    end
    provmatrix=provmatrix./size(genlist,2);
    exchange=(provmatrix+provmatrix')/2;
    dlmwrite([folder1,'\',outputname],exchange,'\t');
    end