function CNull=Temp_Null_Model(C,type,varargin)
% The following MATLAB code is tested on MATLAB R2016a.
% TEMP_NULL_MODEL generates permuted-time models with/without preservance of
% lifetimes.
%
% C - Contact sequence in matrix form. Let nc be the number of contacts. C
% is a NC-by-3 matrix. Each row corresponds to a contact. In particular,
% the first and second entry denotes the two interacted nodes, and the
% thrid entry gives the timestamp of that contact.
%
% type - {'PT'|'PTN'|'PTE'}. We implemented three versions of permuted-time
% reference models. 'PT' is the most common version where timestamps are
% simply reshuffled. In 'PTN', permutaions that change lifetime of any
% node are rejected. Hence, node lifetimes are preserved in 'PTN'. IN
% 'PTE', permutations that change lifetime of any edges are rejected. Both
% node lifetimes and edge lifetimes are retained in 'PTE'.
%
% varargin = {Imax} - The maximum number of permutation attempts for each permutation.
% When the number of attempts is larger than Imax, we move to the next
% permutation. By default, Imax = 100.
%
% CNull - Ramdomzied contact sequence.
%
% Source: M. Li, V.D. Rao, T. Gernat, H. Dankowicz (2018). 
% Lifetime-preserving reference models for characterizing spreading
% dynamics on temporal networks. Scientific Reports 8:709.

%% perform randonization
nc=size(C,1);
CNull=C;
switch type
    case 'PT'
        shuffledinces=randperm(nc);
        CNull(:,3)=C(shuffledinces,3);
    case 'PTE'    
        if ~isempty(varargin) && isa(varargin{1}, 'numeric')
            Imax = varargin{1};
        else
            Imax = 100;
        end        
        % compute birth and death times of edges
        ActInt=Cal_Edge_Life(C);
        ned=size(ActInt,1);
        CNull=C;
        for kk=1:ned
            ii=ActInt(kk,1);
            jj=ActInt(kk,2);
            %% find all contacts for ii-jj or jj-ii
            idx=find(CNull(:,1)==ii);
            id_ij=idx(find(CNull(idx,2)==jj)); %ii-jj
            idx=find(CNull(:,2)==ii);
            id_ji=idx(find(CNull(idx,1)==jj)); %jj-ii
            id=union(id_ij,id_ji);
            %% extract contact times
            tid=CNull(id,3);
            [tid, idd]=sort(tid);
            id=id(idd); % arrange based on time order
            if numel(tid)>=3
                CNull(id(1),3)=tid(1);                                              % preserve the first contact time
                CNull(id(end),3)=tid(end);                                          % preserve the last contact time
                tokid=intersect(find(CNull(:,3)<tid(end)),find(CNull(:,3)>tid(1))); % extract all contact within (tid(1),tid(f))
                for k=2:numel(id)-1
                    nit=0;
                    while 1
                        nit=nit+1;
                        if nit>Imax
                            break;
                        end
                        pid=tokid(randi(numel(tokid),1));                           % generate index for the other contact to permute
                        i1=CNull(pid,1);
                        j1=CNull(pid,2);
                        % find the birth and death times for i1-j1 or j1-i1 edge
                        idx=find(ActInt(:,1)==i1);
                        id_ij=idx(find(ActInt(idx,2)==j1)); %i1-j1
                        idx=find(ActInt(:,2)==i1);
                        id_ji=idx(find(ActInt(idx,1)==j1)); %j1-i1
                        idd=union(id_ij,id_ji);
                        assert(numel(idd)<=1,'multiple edges occured');
                        % find all contacts for i1-j1 or j1-i1 edge
                        idx=find(CNull(:,1)==i1);
                        id_ij=idx(find(CNull(idx,2)==j1)); %ii-jj
                        idx=find(CNull(:,2)==i1);
                        id_ji=idx(find(CNull(idx,1)==j1)); %jj-ii
                        iddd=union(id_ij,id_ji);                                               
                        if CNull(id(k),3)>ActInt(idd,3) && CNull(id(k),3)<ActInt(idd,4)
                            if CNull(pid,3)~=ActInt(idd,3) && CNull(pid,3)~=ActInt(idd,4) 
                                if isempty(find(CNull(id,3)==CNull(pid,3))) && isempty(find(CNull(id(k),3)==CNull(iddd,3))) 
                                    ss=CNull(id(k),3);
                                    CNull(id(k),3)=CNull(pid,3);
                                    CNull(pid,3)=ss;
                                    break;
                                end
                            
                            end
                        end
                    end
                end
            end

        end
        
                 
    case 'PTN'
        if ~isempty(varargin) && isa(varargin{1}, 'numeric')
            Imax = varargin{1};
        else
            Imax = 100;
        end
        % compute birth and death times of nodes
        Life  = Cal_Node_Life(C);
        for k=1:nc
            i=CNull(k,1);
            j=CNull(k,2);
            t=CNull(k,3);
            % compute the intersection of lifetimes of nodes i and j
            idx=find(Life(:,1)==i);
            tai=Life(idx,2);
            tbi=Life(idx,3);
            if Life(idx,4)<3   % node i has only two contacts
                continue
            end
            idx=find(Life(:,1)==j);
            taj=Life(idx,2);
            tbj=Life(idx,3);
            if Life(idx,4)<3   % node j has only two contacts
                continue
            end
            ta=max(tai,taj);
            tb=min(tbi,tbj);
            if t==ta || t==tb  % preserve lifetimes
                continue
            end         
            % find permuted contact randomly
            nit=0;
            while 1
                nit=nit+1;
                if nit>Imax   % upper limit for finding
                    break;
                end
                pset=setdiff(1:nc,k);
                pid=randi(nc-1,1);
                pid=pset(pid);
                ip=CNull(pid,1);
                jp=CNull(pid,2);
                tp=CNull(pid,3);
                if tp>ta && tp<tb
                    % compute the intersection of lifetimes of nodes ip and
                    % jp
                    idx=find(Life(:,1)==ip);
                    tapi=Life(idx,2);
                    tbpi=Life(idx,3);
                    idx=find(Life(:,1)==jp);
                    tapj=Life(idx,2);
                    tbpj=Life(idx,3); 
                    tap=max(tapi,tapj);
                    tbp=min(tbpi,tbpj);   
                    if t>tap && t<tbp
                            if tp~=tap && tp~=tbp  % preserve lifetimes
                                % find all event for i-j edge
                                idx=find(CNull(:,1)==i);
                                id_ij=idx(find(CNull(idx,2)==j)); %i-j
                                idx=find(CNull(:,2)==i);
                                id_ji=idx(find(CNull(idx,1)==j)); %j-i
                                idd=union(id_ij,id_ji); 
                                % find all envent for ip-jp edge
                                idx=find(CNull(:,1)==ip);
                                id_ij=idx(find(CNull(idx,2)==jp)); %ip-jp
                                idx=find(CNull(:,2)==ip);
                                id_ji=idx(find(CNull(idx,1)==jp)); %jp-ip
                                iddd=union(id_ij,id_ji);   
                                % an edge at a timestamp has most one contact
                                if isempty(find(CNull(idd,3)==tp)) && isempty(find(CNull(iddd,3)==t))                                
                                    CNull(k,3)=tp;
                                    CNull(pid,3)=t;
                                    break;
                                end
                            end
                    end
                end
           end
        end      
            
    otherwise
        error('Only PT, PTN and PTE reference models are available.');
end



end







%% ---------------------------------------------------------
function y=Cal_Edge_Life(C)
% Calculate the birth and death times of each edge
% INPUT: C - contact sequence
% OUTPUT: y - lifetime lists of edges. Each row [i j t_0 t_f] corresponds
% to the lifetimes of an edge. i and j denotes the two linked nodes of the
% edge. t_0 and t_f represents the birth and death times of the edge.


y=[];
A=C;
while (~isempty(A))
    alldofs=1:size(A,1);
    ii=A(1,1);jj=A(1,2);
    % find all contacts for ii-jj or jj-ii
    idx=find(A(:,1)==ii);
    id_ij=idx(find(A(idx,2)==jj)); %ii-jj
    idx=find(A(:,2)==ii);
    id_ji=idx(find(A(idx,1)==jj)); %jj-ii
    id=union(id_ij,id_ji);
    % find birth and death times for ii-jj edge
    t_0=min(A(id,3));
    t_f=max(A(id,3));
    y=[y;ii jj t_0 t_f];
    freedofs=setdiff(alldofs,id);
    A=A(freedofs,:);
end

end




%% ------------------------------------------------------------
function y=Cal_Node_Life(C)
% Calculate the birth and death times of each node
% INPUT: C - contact sequence
% OUTPUT: y - lifetime lists of nodes. Each row [i t_0 t_f nc] corresponds
% to the lifetimes of a node with index i. t_0 and t_f represents the 
% birth and death times of the edge. As a supplementary, the number of
% contacts of node i, i.e., nc, is provided as well.

y=[];
node_ids=unique([C(:,1);C(:,2)]);
n=length(node_ids);
for k=1:n
    ii=node_ids(k);
    % find all contacts for ii- or -ii
    id_out=find(C(:,1)==ii); %ii-
    id_in=find(C(:,2)==ii); %-ii
    id=union(id_in,id_out);
    % extrach birth and death times for node ii
    t_0=min(C(id,3));
    t_f=max(C(id,3));
    y=[y;ii t_0 t_f numel(id)];
end


end