 % setup of dns data
 % uip set of unique IP address indexes
 % ipset{i} contains the indexes of the domain-names served by the i-th IP addres of uip
 % pipset(i) probability of finding i-th address of uip in a trace
 % intrapipset{i}[j] probability of having the domain-name j-th given the use of i-th IP address  
 
 load hip.txt
 uip=unique(hip(:,2));
 ipset=cell(length(uip),1);
 for i=1:length(uip)
  ipset(i)={find(hip(:,2)==uip(i))};
 end
 
 fprintf('Number of unique public IP addresses %d \n',length(uip));
 
 alpha=1;
 fprintf('Evaluation of occurrence probability of public IP addresses and domain-names with zipf alpha %d...\n',alpha);
 n=[1:length(hip)];
 z=sum(1./(n.^alpha));
 for i=1:length(uip)
  v=ipset{i};
  pipset(i)=0;
  intrav=[];
  for j=1:length(v)
      p=((1/(v(j)^alpha))/z);
      pipset(i)=pipset(i)+p;
      intrav=[intrav p];
  end
  intrapipset{i}=intrav/pipset(i);
end
 
 
% 
% % the stat_e.txt file should be a copy of the stat-flow-* file wrote by the python modules
% 

sfile=fopen('stat_r.txt','w');
s = RandStream('mt19937ar','Seed', 5811);
RandStream.setDefaultStream(s)

load stat_e.txt
l=length(stat_e(:,1));


% fix the end of routes not ended in the trace

last_full_entry=find(stat_e(:,4)==-1,1)-1;
last_second=max(stat_e(:,2));
last_usecond=max(stat_e(:,3));
stat_e(last_full_entry+1:end,4)=last_second*ones(l-last_full_entry,1);
stat_e(last_full_entry+1:end,5)=last_usecond*ones(l-last_full_entry,1);


stat_e=[stat_e (stat_e(:,2)*1e6+stat_e(:,3)) ((stat_e(:,4)-stat_e(:,2))*1e6+(stat_e(:,5)-stat_e(:,3)))  (stat_e(:,4)*1e6+stat_e(:,5))];


% evaluation of number of active flows 
arrival_events_no=[ones(l,1) stat_e(:,1) stat_e(:,6) stat_e(:,7) (1:length(stat_e))']; %not time ordered
arrival_events=sortrows(arrival_events_no,3); % time ordered
departure_events_no=[-1*ones(l,1) stat_e(:,1) stat_e(:,8) stat_e(:,7) (1:length(stat_e))']; %not time ordered
ev=[arrival_events_no;departure_events_no];
ev=sortrows(ev,3);
ao=cumsum(ev(:,1));
max_index=find(ev(:,3)>ev(end,3)-5e6,1); % remove the last five seconds from the analysis
min_index=find(ev(:,3)>ev(1,3)+5e6,1); % remove the first five seconds from the analysis
ro=mean(ao(min_index:max_index));
max_ao=max(ao(min_index:max_index));
plot(ao)
hold on
fprintf('Average number of active flows of the trace  %d \n',ro);

max_index=find(arrival_events(:,3)>arrival_events(end,3)-5e6,1);
min_index=find(arrival_events(:,3)>arrival_events(1,3)+5e6,1);
fprintf('Average active-flow inter-time %d \n',mean(diff(arrival_events(min_index:max_index,3))));


% unique set of IP addresses of the trace
uip_trace=unique(arrival_events(:,2));
fprintf('Number of unique anonymous IP addresses of the trace %d \n',length(uip_trace));

% evaluation of recurrence probability of each unique IP ananymous address of the
% trace
fprintf('Evaluation of occurrence probability of anonymous IP addresses ... \n');
n_flow_uip_trace = histc(arrival_events(:,2),uip_trace);
n_flow_uip_trace = n_flow_uip_trace / sum(n_flow_uip_trace);
v_uip_trace=[uip_trace n_flow_uip_trace [1:length(uip_trace)]'];
v_uip_trace=sortrows(v_uip_trace,2);

fprintf('Simulation start ... \n');
pipset_temp=[pipset' [1:length(pipset)]'];
% (dottedQuadToNum(key),first_arrival_s, first_arrival_us, last_arrival_s, last_arrival_us

fprintf('Random extraction of %d public IP addresses based on their occurrence probability ...\n',length(uip_trace));
rv=rand(length(pipset),1);
cs=cumsum(pipset_temp(:,1)/sum(pipset_temp(:,1)));
r=rv(1);
index_uip_dns_temp=dsearchn(cs,r);
if(cs(index_uip_dns_temp)<r) 
   index_uip_dns_temp=index_uip_dns_temp+1;
end
index_uip_dns=pipset_temp(index_uip_dns_temp,2);
v_uip_dns=zeros(length(uip_trace),3);
v_uip_dns(1,:)=[uip(index_uip_dns) pipset(index_uip_dns)' index_uip_dns];
pipset_temp(index_uip_dns_temp,:)=[];

for i=2:length(uip_trace)
   cs=cumsum(pipset_temp(:,1)/sum(pipset_temp(:,1)));
   r=rv(i);
   index_uip_dns_temp=dsearchn(cs,r);
   if(cs(index_uip_dns_temp)<r) 
       index_uip_dns_temp=index_uip_dns_temp+1;
   end
   index_uip_dns=pipset_temp(index_uip_dns_temp,2);
   v_uip_dns(i,:)=[uip(index_uip_dns) pipset(index_uip_dns)' index_uip_dns];
   pipset_temp(index_uip_dns_temp,:)=[];
   
   if (mod(i,1000)==0)
       i
   end
end
v_uip_dns=sortrows(v_uip_dns,2);
if (length(unique(v_uip_dns(:,3)))<length(uip_trace))
     fprintf('Extraction fail: number of extracted public addresses lower that anonymized IP addresses');
     return;
end

fprintf('Run time association of anonymous flows to domain-names ... \n');

fib_size=10000; %% this should be set so large as the FIB was infinite
s=ones(fib_size,5)*-1;   % vector containig the id of the route contained in fib entry, the next free time of fib entry and the start time 
arrived=0;
min_now=arrival_events(1,3);
arrival_events(:,3)=arrival_events(:,3)-min_now;
last_sample_time=0;
last_sample_time_stamp=0;
ito=1e6;
for i=1:length(arrival_events)   %max_index

    % Assigning a domain-name to a flow  
    ip_trace=arrival_events(i,2);
    ip_trace_index=find(v_uip_trace(:,1)==ip_trace);
    ip_dns_index=v_uip_dns(ip_trace_index,3);
    cs=cumsum(intrapipset{ip_dns_index});
    r=rand(1);
    index_intra=dsearchn(cs',r);
    if(cs(index_intra)<r) 
           index_intra=index_intra+1;
    end
    v=ipset{ip_dns_index};
    domain_name_index=v(index_intra);
    % remapping complete


    now=arrival_events(i,3);
    si=find(s(:,1)==domain_name_index,1);
    if (numel(si)>0)
      % ICN route already contained in the FIB
      % refresh expiry time s(si,2)
      s(si,3)=now;
      if (now+arrival_events(i,4)>s(si,2))
        s(si,2)=now+arrival_events(i,4);
        s(si,5)=arrival_events(i,5);
      end
    else
        % ICN route not contained in the FIB 
        [min_value,min_i]=min(s(:,2));  % select a free entry (s(j,2)==-1)
        if (min_value==-1)
            s(min_i,1)=domain_name_index;
            s(min_i,2)=now+arrival_events(i,4);
            s(min_i,3)=now;
            s(min_i,4)=arrival_events(i,5);
            s(min_i,5)=arrival_events(i,5);
        else
            fprintf('Error critical route blocking, FIB size not enough !');
            return;
        end            
    end
    if (now-last_sample_time>100e3) % sampling period about 100e3 ms
        
        for k=1:fib_size
            if (s(k,1)>-1 && s(k,2)+ito<now)
                % (dottedQuadToNum(key),first_arrival_s, first_arrival_us, last_arrival_s, last_arrival_us
                asec=stat_e(s(k,4),2);%floor(s(k,3)/1e6);
                ausec=stat_e(s(k,4),3);%(s(k,3)-asec*1e6);
                dsec=stat_e(s(k,5),4);%floor(s(k,2)/1e6);
                dusec=stat_e(s(k,5),5);%(s(k,2)-dsec*1e6);
                fprintf(sfile,'%d %d %d %d %d \n', s(k,1),asec,ausec,dsec,dusec);
                s(k,:)=[-1 -1 -1 -1 -1];
            end
        end
        last_sample_time=now;
        if (now-last_sample_time_stamp>1e6)
            nactiveroutes=length(find(s(:,2)>=now));
            fprintf('\n Estimted number of ICN active-routes %d at %d\n',nactiveroutes,now/1e6);  
            last_sample_time_stamp=now;
        end
    end
end
% finalize remaining active routes
for k=1:fib_size
    if (s(k,1)>-1)
        % (dottedQuadToNum(key),first_arrival_s, first_arrival_us, last_arrival_s, last_arrival_us
        asec=stat_e(s(k,4),2);%floor(s(k,3)/1e6);
        ausec=stat_e(s(k,4),3);%(s(k,3)-asec*1e6);
        dsec=stat_e(s(k,5),4);%floor(s(k,2)/1e6);
        dusec=stat_e(s(k,5),5);%(s(k,2)-dsec*1e6);
        fprintf(sfile,'%d %d %d %d %d \n', s(k,1),asec,ausec,dsec,dusec);
        s(k,:)=[-1 -1 -1 -1 -1];
    end
end
fclose(sfile);

