cr_bf.m

best routing protocol

function  [routing_dist, routing_load, storage_load, dht_packets, cluster] = ...
    cr_bf(neighbor_dist,pre_landmarks, ...
		   test_pairs, failure_nodes, reactive)
%Refer to "Compact routing schemes" by Thorup & Zwick for the
%algorithm 
%neighbor_dist: distances are arbitrary ones, not necessarily
%hops. inf indicates non-connectivity.
% 
% pre_landmarks: pre-set landmark sets. 
% if not given, used the proposed randomized algorithm in the paper
% if length(pre_landmarks)==1, it's used as the number of landmarks

debug = 1;
debug_route = 0;
if ~issparse(neighbor_dist)
  neighbor_dist = dist_full2sparse(neighbor_dist);
end
if nargin<5
  reactive = 1; %default CR is proactive, i.e. reactive = 0;
end

n = length(neighbor_dist);
%fprintf('calculating best path...');
%[delta optimal_nexthop] = floyd_warshall(neighbor_dist);
%[delta optimal_nexthop] = allpair_dijkstra(neighbor_dist);
%fprintf('done\n');

%%%select landmarks
W = 1:n;
if nargin<2
  pre_landmarks = []; 
end

if nargin<3
  tmp1 = repmat(1:n,n,1);
  tmp2 = tmp1';
  rest_pairs = [tmp1(:) tmp2(:)];
end

if nargin<4
  failure_nodes = [];
end

fprintf('CR: reactive mode:%d\n',reactive);

out_degree = full(sum(neighbor_dist>0,2));
%C = zeros(n,n); %C(w,v)==1 means v \in C(w)
%routing_table(w,v)=nexthop, if v \in C(w)
routing_table = sparse([]);
routing_table(n,n) = 0;
direct_dist = routing_table; %should be sync with routing_table:
                             %direct_dist(w,v)=shortest_dist bwteen w&v

%nexthop2bcn(i,b) = nexthop(i,b) where b is the relative beacon id
deltaA = zeros(1,n); % deltaA(v) == \delta(A,v)
%cent(u), the closest landmark to each node, relative index to A

A = [];
s = sqrt(n*log(n)/log(2));
%s = sqrt(n);%!!!!
while length(W)>0
  if length(pre_landmarks)>=1%!!!!
    A = pre_landmarks;
  elseif length(pre_landmarks)==1
    % pre_landmarks is defined as the # of random landmarks
    tmp = randperm(n);
    A = tmp(1:pre_landmarks);
  else
    A = union(A, sample(W,s));
  end
  A = sort(A);
  ordinary_nodes = setdiff(1:n, A);
  isbeacon = zeros(n,1);
  isbeacon(A) = 1:length(A);
  %every landmark floods the network
  lm_broadcast_snd = length(A)*n;
  lm_broadcast_rcv = length(A)*sum(out_degree);
  
  if debug
    fprintf('# broadcast packets sent: %d\n', lm_broadcast_snd);
    fprintf('# broadcast packets received: %d\n', lm_broadcast_rcv);
  end

  [beacon_vectors, nexthop2bcn] = allpair_dijkstra(neighbor_dist, A);
  [deltaA, cent] = min(beacon_vectors,[],2);
  %plot_cumu_matrix(deltaA);%!!!!
  if debug
    fprintf('mean and max deltaA:%.1f, %.1f\n', mean(deltaA), ...
	    max(deltaA));
  end
  if debug & n<=50
    deltaA'
  end
  
  %%%now get treeLabel of each node
  fprintf('CR: calculating tree labels...\n');
  tlabel = cell(n,1);
  for i=ordinary_nodes
    bcn = cent(i);
    if bcn<=0
      error('check here');
    end
    
    if deltaA(i)==Inf
      tlabel{i}(1) = -1;
    else
      curr = i;
      for count = deltaA(i)-1:-1:1
        curr = nexthop2bcn(curr, bcn);
        tlabel{i}(count) = curr;
      end
    end
  end
  %%%
  
  %now sort out the clusters for ordinary nodes
  %%C = delta < repmat(deltaA,n,1);
  %%C = delta <= repmat(deltaA,n,1); %!!!change from the algorithm
  %%C(:,A) = zeros(n, length(A));
  %%C(A,:) = zeros(length(A),n);
  if reactive
    fprintf('do NOT do local broadcasting\n')
  else
    fprintf('simulating local broadcasting...\n');
    % every ordinary node i floods the network in the radix of deltaA(i) to let others find itself in their clusters
    od_broadcast_snd = 0;
    od_broadcast_rcv = 0;
    
    count = 0; total_c = length(ordinary_nodes); old_c = 0;
    start_t = clock;old_t = clock;
  
    for i=ordinary_nodes
      new_t = clock;
      if etime(new_t, old_t)>10
	fprintf('progress %.1f%%. ETA %.1f secs\n', count*100/total_c, ...
		etime(new_t,old_t)/(count-old_c)*(total_c-count));
	old_t = new_t;
	old_c = count;
      end
      
      [b_snd, b_rcv, rt_i,label(i),rt_dist ] = BFS(i,deltaA(i), neighbor_dist, ...
					   1, A(cent(i)) );
      routing_table(:,i) = rt_i';
      direct_dist(:,i) = rt_dist';
      
      routing_table(n,n) = 0;
      od_broadcast_snd = od_broadcast_snd + b_snd;
      od_broadcast_rcv = od_broadcast_rcv + b_rcv;
      count = count+1;
    end
  
    if debug
      fprintf('# broadcast packets sent due to local clustering: %d\n', od_broadcast_snd);
      fprintf('# broadcast packets received due to local clustering: %d\n', od_broadcast_rcv);
    end
  end
  %%%%
  if length(pre_landmarks)>0
    W = [];
  else
    W = find( sum(C,2) > 4*n/s);
  end
  
end
if debug
  fprintf('%d selected landmarks.\n',length(A));
  if n<100
    A
  end
end
%% end of selecting landmarks.
%% other useful variables: cent(dest), C(src,dest)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% use bloom filter to store routing information

%% merge routing_table, next2bcn and tlabel 

% rtb{i} is the routing table at node i
% each entry has two elements, dest and nexthop 
rtb = cell(n,1);

% from routing_table
for i=1:n
  v = find(routing_table(i,:)>0)';
  rtb{i} = [rtb{i}; [v routing_table(i,v)']];
end

% from nexthop2bcn
for i=1:n
  b = find(nexthop2bcn(i,:)>0)';
  rtb{i} = [rtb{i}; [A(b)' nexthop2bcn(i,b)']];
end

% from tlabel
for i=1:n
  if ~isempty(tlabel{i}) && tlabel{i}(1) > 0
    b = A(cent(i));	% landmark closest to i
    rtb{b} = [rtb{b}; [i tlabel{i}(1)]];
  end
end

%% build routing bloom filters

% rbf{i} is the bloom filter at node i
rbf = cell(n,1);

% adapted from Universal Class Hash Functions
ratio = 15;	% ratio of the bit filter size to the numbet of keys
hk = 8;		% number of independent hash functions
%p = 10007;	% smallest prime number greater than 10000
p = [10007 10009 10037 10039 10061 10067 10069 10079];
%coe1 = floor(rand([1 hk])*(p-2))+1;	0<coe1<p
%coe2 = floor(rand([1 hk])*(p-1));	0<=coe2<p
coe1 = [2092 3801 7838 6812 4614 5682 7947 593];
coe2 = [6032 502 4156 3051 8748 150 7684 9714];

for i=1:n
  %rtb{i}
  bitlength = size(rtb{i},1)*ratio;
  rbf{i} = zeros(bitlength,1);
  for j=1:hk
    bitpos = mod(mod(coe1(j)*rtb{i}(:,1) + coe2(j)*rtb{i}(:,2),p(j)),bitlength);
    rbf{i}(bitpos+1) = 1;
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%use DHT to store node coordinates, see the overhead
%storage overhead per beacon, + packets sent
fprintf('CR: putting vectors to beacons as a DHT\n');
beacons = A;
ordinary_nodes = setdiff(1:n, beacons);
storage_load = zeros(length(beacons),1);
dht_packets = 0;
for i=ordinary_nodes
  b_k = consistent_hash(i,beacons, n); %b_k \in 1.. length(beacons)
  storage_load(b_k) = storage_load(b_k)+1;
  tmp_dist = 0;
  [f_path, f_dist, routing_table, direct_dist] = cr_routing(i, ...
					       A(b_k), cent, A, debug_route, beacon_vectors, neighbor_dist, ...
					       routing_table, nexthop2bcn, isbeacon, [],reactive, ...
					       deltaA, direct_dist, ...
						  failure_nodes);
  
  [b_path, b_dist, routing_table, direct_dist] = cr_routing(A(b_k), ...
					       i, cent, A, debug_route, beacon_vectors, neighbor_dist, ...
					       routing_table, nexthop2bcn, isbeacon, tlabel{i},reactive, ...
					       deltaA, direct_dist, ...
						  failure_nodes);
  dht_packets = dht_packets + f_dist + b_dist;
end
if length(beacons)<10
  storage_load
end

fprintf('CR: dht_packets:%d; mean storage load:%.2f, max:%.2f\n', ...
	 dht_packets, mean(storage_load), max(storage_load));
%%%%


%% routing algorithm
if 0 %old.
  routing_nexthop = optimal_nexthop(:, A(cent)); %route to the
						 %nearest landmarks