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(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)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 && tptap && t