diff --git a/src/DHCMall.m b/src/DHCMall.m
new file mode 100644
index 0000000..a12f834
--- /dev/null
+++ b/src/DHCMall.m
@@ -0,0 +1,37 @@
+function [f,g]=DHCMall(G,N)
+n=length(G);
+%
+measureA=cell(N,1); measureB=measureA;
+%[C1,C2,C]=clust_coeff(G);
+measureA{1}=G; measureB{1}=G;
+A=G; B=G;
+X=zeros(n,N);
+%Standard clustering:
+for i=1:n
+    vi=A(i,:);
+    vfind=find(vi);
+    lv=length(vfind);
+    Wi=B(vfind,vfind);
+    X(i,1)=sum(sum(Wi))/lv/(lv-1);%length(find(Wi))/lv/(lv-1);
+end
+%N-clustering:
+for m=2:N
+    nextA=A*G;
+    nextA(speye(n)==1)=0;
+    nextA(B==1)=0;
+    nextA(nextA>1)=1;
+    nextB=B+nextA;
+    nextB(speye(n)==1)=0;
+    nextB(nextB>1)=1;
+    measureA{m}=nextA; measureB{m}=nextB;
+    A=nextA; B=nextB;
+    for i=1:n
+        vi=B(i,:);%Ax if exactly N steps from i, Bx if up to N steps
+        vfind=find(vi);
+        lv=length(vfind);
+        Wi=B(vfind,vfind);%B1 for neighbours, Bx for up to N steps between, BNm1 for N-1
+        X(i,m)=sum(sum(Wi))/lv/(lv-1);%length(find(Wi))/lv/(lv-1);
+    end
+end
+X=X';
+f=[nanmean(X,2),nanvar(X,0,2),min(X,[],2),prctile(X,[5,25,50,75,95],2),max(X,[],2)];
\ No newline at end of file
diff --git a/src/DHCMsampling.m b/src/DHCMsampling.m
new file mode 100644
index 0000000..182677e
--- /dev/null
+++ b/src/DHCMsampling.m
@@ -0,0 +1,107 @@
+function f=DHCMsampling(G,mmax,sampleNumber)%Or input v
+%G must be simple: undirected with no self loops
+n=length(G);
+v=randsample(n,sampleNumber,false);%v is indices
+lv=length(v);
+clusterVec1=zeros(lv,mmax);
+
+imvw1=v';%node
+imvw2=ones(lv,1);%info order
+Cv=cell(n,1);%{};%Cell array of data for each node
+lengthC=lv;
+
+%Standard clustering:
+for i=1:lv
+    %Find neighbours of node i:
+    vx=v(i);
+    v1=mNbr(G,vx,1,vx);%Indices of neighbours
+    %Update stored data:
+    Cv{i}=v1;
+    %
+    Gred=G(v1,v1);%Adj mat of neighbourhood
+    lv1=length(v1);%Number of neighbours
+    C1=sum(sum(Gred))/lv1/(lv1-1);%How many links exist/how many could exist
+    clusterVec1(i,1)=C1;
+end
+%m-clustering, m>1:
+if mmax>1
+for m=2:mmax
+    for i=1:lv
+        %m-neighbours of v(i):
+        vi=v(i);
+        vx=Cv{i};%m-1 OR m (if already xcalculated in loop) nbrs of v(i)
+        imvw2i=imvw2(i);
+        if imvw2i<m%i.e.==m-1
+            vm=mNbr(G,vx,m-imvw2i,vi);%Neighbours of i
+            %vm=vm(vm~=vi) here?
+            imvw2(i)=m;
+            Cv{i}=vm;
+        else
+            vm=vx;
+        end
+        vmx=vm;%(vm~=vi);%Exclude node i
+        lvm=length(vmx);
+        %m-clustering:
+        %
+        links1=0;
+        %m-neighbours of each m-neighbour:
+        for j=1:lvm%parfor
+            vmj=vmx(j);%vmOrder(j);
+            %if vmj~=vi
+            %Faster to order or to find all links?
+            if ismember(vmj,imvw1)%If have already calculated some nbhds
+                index=find(imvw1==vmj);%How many
+                thisFar=imvw2(index);
+                vmjSoFar=Cv{index};
+                if thisFar<m
+                    mNj=mNbr(G,vmjSoFar,m-thisFar,vmj);
+                    imvw2(index)=m;
+                    Cv{index}=mNj;
+                else
+                    mNj=Cv{index}; %mNjOnly=Cw{index};
+                end
+            else
+                mNj=mNbr(G,vmj,m,vmj);
+                lengthC=lengthC+1;
+                index=lengthC;
+                imvw1(index)=vmj;
+                imvw2(index)=m;
+                Cv{index}=mNj;
+            end
+            mNjx=mNj;%(mNj~=vmj);%m-neighbours of node j, not including j
+            links1=links1+length(intersect(vmx,mNjx));
+        end
+        clusterVec1(i,m)=links1/lvm/(lvm-1);
+    end
+end
+end
+%****Save output****
+%f=clusterVec1;
+X=clusterVec1';
+f=[nanmean(X,2),nanvar(X,0,2),min(X,[],2),prctile(X,[5,25,50,75,95],2),max(X,[],2)];
+end
+
+function f=mNbr(G,vx1,mMore,vi)
+    vx=vx1;%Neighbours so far (to order m-mmore)
+    G1=G(vx,:);%Next neighbours (1)
+    sumG1=sum(G1,1);%Vector next neighbours (>=1, accounting for overlap) (2)
+    v1=find(sumG1>0);%Indices of next neighbours (3)
+    
+    v1=union(vx,v1);%up-to-m-neighbours of i (4)
+    
+    v1=v1(v1~=vi);%Exclude original node (5)
+    vx=v1;%In case want wx etc.
+    if mMore>1
+        k=2;
+        while k<=mMore
+            Gm=G(vx,:);%Repeat (1)
+            sumGm=sum(Gm,1);%(2)
+            vm=find(sumGm>0);%(3)
+            vm=union(vx,vm);%up-to-m-neighbours of i (4)
+            vm=vm(vm~=vi);%(5)
+            vx=vm;
+            k=k+1;
+        end
+    end
+    f=vx;
+end
\ No newline at end of file
diff --git a/src/DHCMsamplingSpeed.m b/src/DHCMsamplingSpeed.m
new file mode 100644
index 0000000..5de0539
--- /dev/null
+++ b/src/DHCMsamplingSpeed.m
@@ -0,0 +1,114 @@
+function f=DHCMsamplingSpeed(G,mmax,sampleNumber)%Or input v
+%G must be simple: undirected with no self loops
+n=length(G);
+v=randsample(n,sampleNumber);%v is indices
+%%
+lv=length(v);
+clusterVec1=zeros(lv,mmax);
+
+%Gnbrs=cell{n,1};
+[aye,jay]=find(G);
+GG=accumarray(aye,jay,[size(G,1),1],@(x) {sort(x).'});
+
+%%
+imvw1=v';%node
+imvw2=ones(lv,1);%info order
+Cv=cell(n,1);%{};%Cell array of data for each node
+lengthC=lv;
+
+%Standard clustering:
+for i=1:lv
+    %Find neighbours of node i:
+    vx=v(i);
+    v1=mNbr(GG,vx,1,vx);%Indices of neighbours
+    %Update stored data:
+    Cv{i}=v1;
+    %
+    Gred=G(v1,v1);%Adj mat of neighbourhood
+    lv1=length(v1);%Number of neighbours
+    C1=sum(sum(Gred))/lv1/(lv1-1);%How many links exist/how many could exist
+    clusterVec1(i,1)=C1;
+end
+%m-clustering, m>1:
+if mmax>1
+for m=2:mmax
+    for i=1:lv
+        %m-neighbours of v(i):
+        vi=v(i);
+        vx=Cv{i};%m-1 OR m (if already xcalculated in loop) nbrs of v(i)
+        imvw2i=imvw2(i);
+        if imvw2i<m%i.e.==m-1
+            vm=mNbr(GG,vx,m-imvw2i,vi);%m-neighbours of i
+            %vm=vm(vm~=vi) here?
+            imvw2(i)=m;
+            Cv{i}=vm;
+        else
+            vm=vx;
+        end
+        vmx=vm;%(vm~=vi);%Exclude node i
+        lvm=length(vmx);
+        %m-clustering:
+        %
+        links1=0;
+        %m-neighbours of each m-neighbour:
+        for j=1:lvm%parfor
+            vmj=vmx(j);%vmOrder(j);
+            %if vmj~=vi
+            %Faster to order or to find all links?
+            if ismember(vmj,imvw1)%If have already calculated some nbhds
+                index=find(imvw1==vmj);%How many
+                thisFar=imvw2(index);
+                vmjSoFar=Cv{index};
+                if thisFar<m
+                    mNj=mNbr(GG,vmjSoFar,m-thisFar,vmj);
+                    imvw2(index)=m;
+                    Cv{index}=mNj;
+                else
+                    mNj=Cv{index}; %mNjOnly=Cw{index};
+                end
+            else
+                mNj=mNbr(GG,vmj,m,vmj);
+                lengthC=lengthC+1;
+                index=lengthC;
+                imvw1(index)=vmj;
+                imvw2(index)=m;
+                Cv{index}=mNj;
+            end
+            %Other method: find subgraph
+            %mNj=mNj(~vi); mNjOnly=mNjOnly(~vi);
+            mNjx=mNj;%(mNj~=vmj);%m-neighbours of node j, not including j
+            links1=links1+length(intersect(vmx,mNjx));
+            %end
+        end
+        clusterVec1(i,m)=links1/lvm/(lvm-1);
+    end
+end
+%****Save output****
+%f=clusterVec1;
+end
+X=clusterVec1';
+f=[nanmean(X,2),nanvar(X,0,2),min(X,[],2),prctile(X,[5,25,50,75,95],2),max(X,[],2)];
+end
+
+function f=mNbr(GG,vx1,mMore,vi)
+    vx=vx1;
+    G1cell=GG(vx);
+    v1=unique([G1cell{:}]);%Indices of next neighbours
+    
+    v1=union(vx,v1);%up-to-m-neighbours of i
+    
+    v1=v1(v1~=vi);%Exclude original node
+    vx=v1;%In case want wx etc.
+    if mMore>1
+        k=2;
+        while k<=mMore
+            Gmcell=GG(vx);
+            vm=unique([Gmcell{:}]);
+            vm=union(vx,vm);%up-to-m-neighbours of i
+            vm=vm(vm~=vi);
+            vx=vm;
+            k=k+1;
+        end
+    end
+    f=vx;
+end
\ No newline at end of file
diff --git a/src/RPcellMean.m b/src/RPcellMean.m
new file mode 100644
index 0000000..19d8028
--- /dev/null
+++ b/src/RPcellMean.m
@@ -0,0 +1,13 @@
+function f=RPcellMean(c1,c2,c3)
+%ncells=3;
+%lc=length(c1);
+[a,b]=size(c1{1});
+x=zeros(a,b,ncells);
+cout=cell(lc,1);
+for i=1:lc
+    x(:,:,1)=c1{i};
+    x(:,:,2)=c2{i};
+    x(:,:,3)=c3{i};
+    cout{i}=nanmean(x,3);
+end
+f=cout;
\ No newline at end of file
diff --git a/src/RPcombine.m b/src/RPcombine.m
new file mode 100644
index 0000000..0e0761b
--- /dev/null
+++ b/src/RPcombine.m
@@ -0,0 +1,10 @@
+function f=RPcombine(G)%,fileName)%(B,W)%between(I,J), within(I,J)
+%G=[B(2:end,:);W];
+%lg=length(G);
+lg=max(max(G));
+Gsp=sparse(G(:,1),G(:,2),1,lg,lg);
+Gsp=Gsp+Gsp';
+Gsp(Gsp>1)=1;
+Gsp=Gsp-diag(diag(Gsp));
+%save(fileName,'G')
+f=Gsp;
\ No newline at end of file
diff --git a/src/RPdiffAlphaProcessEpis.m b/src/RPdiffAlphaProcessEpis.m
new file mode 100644
index 0000000..6e0b1fb
--- /dev/null
+++ b/src/RPdiffAlphaProcessEpis.m
@@ -0,0 +1,46 @@
+function [f,g,h]=RPdiffAlphaProcessEpis
+a=(0:1:6);%Alpha values
+n=10000; lruns=10;
+la=length(a);
+P=zeros(la,lruns); Z=P;%5=max poss runs;
+epi=cell(la,1);
+%
+lepi=300;%(Max) length of incidence vector
+epiz=zeros(lepi,lruns);
+[epi{:}]=deal(epiz);
+fpath='..\..\..\Downloads\SRebola\id_spatial_sim-master\scenarios\ebola\output\epiDHh4w30pp2333r2p2_alpha';
+for i=1:la
+    ai=a(i);
+    %events=importdata(strcat('DH_h4w120pp05_alpha',num2str(ai),'_pset_0_Events.out'));%,' ',1,0);
+    events=importdata(fullfile(strcat(fpath,num2str(ai),'_pset_0_Events.out')));
+    %DH1: DH4_30_p3_2_alpha
+    %DH2: DH_h4w100pp08_alpha
+    %DH_h4w120pp05_alpha
+    %
+    run=events.data(:,1);%events(2:end,1);
+    day=events.data(:,2);%events(2:end,2);
+    event=events.data(:,3);%events(3:10:end);%(2:end,3);
+    %index=events(2:end,4);
+    runs=unique(run);
+    lruns=length(runs);
+    for j=1:lruns
+        runj=runs(j);
+        d=day(run==runj)+1;
+        ev=event(run==runj);
+        lev=length(ev); zev=zeros(lev,1);
+        zev(ev==0)=1;
+        %zev=ev; zev(zev>0)=1;
+        x=accumarray(d,zev);
+        [peak,~]=max(x);
+        z=sum(zev);
+        P(i,j)=peak/n; Z(i,j)=z/n;
+        %For full epidemics:
+        epi{i}(1:length(x),j)=x;
+    end
+    zevAll=event; zevAll(zevAll>0)=1; dayAll=day+1;
+    x=accumarray(dayAll,zevAll);
+    epi{i}=x/n/lruns;
+    clear events zevAll dayAll
+end
+f=P; g=Z; h=epi;
+save('allEpis.mat','P','Z','epi')
\ No newline at end of file
diff --git a/src/RPdiffAlphaProcessNets.m b/src/RPdiffAlphaProcessNets.m
new file mode 100644
index 0000000..7006fb5
--- /dev/null
+++ b/src/RPdiffAlphaProcessNets.m
@@ -0,0 +1,45 @@
+function f=RPdiffAlphaProcessNets
+a=(0:6);%Alpha values
+%ai=0;
+%w=[60,100,140,180,220];
+%lw=length(w);
+n=50000;%
+la=length(a);
+reps=1;%Num of nets for each alpha;
+Gcell=cell(la,reps);%XX
+fpath='C:\Users\dhaw\Documents\GitHub\id_spatial_sim\subExp\ebola\network\netDHflath3w200pp04A_alpha';
+for i=1:la
+    ai=a(i);%XX
+    for j=1:reps
+    %
+        fname1=strcat(fpath,num2str(ai),'_arcs.out');
+        fname2=strcat(fpath,num2str(ai),'_nodes.out');
+        arcs=load(fullfile(fname1)); %Households
+        arcs=arcs(:,1);%1st column
+        nodes=load(fullfile(fname2));
+        nodes=nodes(:,5);%5th column
+        %DH1: net_h3w32pp25_alpha
+        %DH2: netDH4_100_p08_2_alpha
+        %DH3: netDH4_120_p05_2_alpha
+        %DH4: netDH4_20_p35_2_alpha (h=3)
+        Gcell{i,j}=RPmakeNets(arcs,nodes,n);
+    end
+    %
+    %{
+    fname1a=strcat(fpath,num2str(ai),'w',num2str(wi),'a_arcs.out');
+    fname2a=strcat(fpath,num2str(ai),'w',num2str(wi),'a_nodes.out');
+    fname1b=strcat(fpath,num2str(ai),'w',num2str(wi),'b_arcs.out');
+    fname2b=strcat(fpath,num2str(ai),'w',num2str(wi),'b_nodes.out');
+    fname1c=strcat(fpath,num2str(ai),'w',num2str(wi),'c_arcs.out');
+    fname2c=strcat(fpath,num2str(ai),'w',num2str(wi),'c_nodes.out');
+    arcs=load(fullfile(fname1a)); arcs=arcs(:,1); nodes=load(fullfile(fname2a)); nodes=nodes(:,5);
+    Gcell{i,1}=RPmakeNets(arcs,nodes,n);
+    arcs=load(fullfile(fname1b)); arcs=arcs(:,1); nodes=load(fullfile(fname2b)); nodes=nodes(:,5);
+    Gcell{i,2}=RPmakeNets(arcs,nodes,n);
+    arcs=load(fullfile(fname1c)); arcs=arcs(:,1); nodes=load(fullfile(fname2c)); nodes=nodes(:,5);
+    Gcell{i,3}=RPmakeNets(arcs,nodes,n);
+    %}
+end
+f=Gcell;
+%save('allNetsDH4.mat','Gcell','-v7.3')%h_w_pw_R0
+save('DH2netsA.mat','Gcell','-v7.3')%h_w_pw_R0
\ No newline at end of file
diff --git a/src/RPmclusPlotsError.m b/src/RPmclusPlotsError.m
new file mode 100644
index 0000000..0a89618
--- /dev/null
+++ b/src/RPmclusPlotsError.m
@@ -0,0 +1,62 @@
+function f=RPmclusPlotsError(X)
+%Input is RPcellMean(cCell)
+%numClus=1;
+clusmax=4;
+la=length(X);
+m1=zeros(clusmax,la);
+M1=zeros(clusmax,la,2);
+for i=1:la
+    Xi=X{i};
+    Xi(:,5)=Xi(:,6)-Xi(:,5);
+    Xi(:,7)=Xi(:,7)-Xi(:,6);
+    M1(:,i,:)=Xi(:,[5,7]);
+    m1(:,i)=Xi(:,6);
+end
+x=repmat((0:6),clusmax,1);
+
+cmap=lines(clusmax);
+%cmap=[.2,.2,.2;1,.2,.2;.2,1,.2;.2,.2,1];
+%cmap=.8*[0,0,0;1,0,0;0,1,0;0,0,1];
+colors=cell(clusmax,1);
+for i=1:clusmax
+    colors{i}=cmap(i,:);
+end
+lineProps.col=colors;
+lineProps.style='-';
+%
+fs=12; lw=2;
+figure
+%lineProps=[];
+%lineProps.col=lines(clusmax);
+hold on
+%
+x1=0:6;
+%
+mseb(x,m1,M1,lineProps,1);
+h1=scatter(x1,m1(1,:),[],cmap(1,:),'o','linewidth',lw);
+h2=scatter(x1,m1(2,:),[],cmap(2,:),'x','linewidth',lw);
+h3=scatter(x1,m1(3,:),[],cmap(3,:),'^','linewidth',lw);
+h4=scatter(x1,m1(4,:),[],cmap(4,:),'s','linewidth',lw);
+hold off
+set(gca,'fontsize',fs)
+axis ([0,6,0,.7])%.35])
+xlabel('Distance power \alpha')
+ylabel(strcat('CC^m (nodes)'))
+legend([h1,h2,h3,h4],'m=1','m=2','m=3','m=4','location','NW')
+%legend([h1,h2],'m=1','m=2','location','NW')
+grid on
+grid minor
+box on
+%}
+%{
+figure
+mseb(x2,h1',H1,lineProps,1);
+set(gca,'fontsize',12)
+axis ([0,6,0,.35])
+xlabel('Distance power \alpha')
+ylabel(strcat('CC^m (households)'))
+legend('m=1','m=2','location','NW')
+grid on
+grid minor
+box on
+%}
\ No newline at end of file