study_cr.m

best routing protocol

function res =  study_cr(n_nodes, density, n_beacons, failure_percent)
% study compact routing, and compare it to optimal routing and BVR
% n_nodes: number of nodes, or the name of .mat file containing the
% whole configuration
% randomly generated topology, density = # of nodes per unit area
%
% usage examples:
% study_cr(1000, 5.12) : simulate 1000 node, 5.12
% density. sqrt(1000) beacons, no failure node

% study_cr('topo-100-5.12-0.1.mat')
% use a pre-generated topology file. 

% study_cr('topo-100-5.12-0.1.mat', 0, 20)
% use a pre-generate topology file, but use different number of
% beacons (20 instead of 10)

% to generate a topology files, first start study_cr() with the
% topology you want, then it's automatically saved in topo-tmp.mat
% rename it to whatever you want.

%%%%%%%%%specify what algorithms you want to simulate here%%%
run_bvr = 1; %simulate bvr or not
run_bvr_edist = 0; % simulate bvr with euclidean distance as dist func
run_optimal = 1;%simulate dijkstra or not
show_all_figures = 0;% plot all the results or not
run_cr_reactive = 0; % simulate the reactive version of CR or not
stretch_percentile = 90;
DEFAULT_NODE_NO = 3200;
%%%%%%end of configuration%%%%%%%%%%%%%%
rand('state',100*(sum(clock)+cputime));
if nargin<1
  n_nodes = DEFAULT_NODE_NO; %this is what is used in the BVR paper
  density = 5.12;
end
topology_fn = '';
if ischar(n_nodes)
  topology_fn = n_nodes;
end
if topology_fn & nargin==1 %load a pre-set configuration
  fprintf('loading everything from %s\n', topology_fn);
  load(topology_fn)
else
  if topology_fn
    %re-specified n_beacons, or failure_percent
    fprintf('loading D & coord from %s\n', topology_fn);
    S = load(topology_fn);
    D = S.D;
    coord = S.coord;    n_nodes = S.n_nodes;
    n_random_pairs = S.n_random_pairs;
    random_pairs = S.random_pairs;
    clear S;
  else
    n_random_pairs = min(ceil(n_nodes*n_nodes/2),DEFAULT_NODE_NO*10);
    if nargin<1
      n_random_pairs = DEFAULT_NODE_NO*10;
      k = 10;
    end
    good_nodes = 1:n_nodes;
    random_pairs_idx = ceil(rand(n_random_pairs, 2)*length(good_nodes));
    random_pairs = reshape(good_nodes(random_pairs_idx(:)), n_random_pairs, ...
			   2);
    clear random_pairs_idx good_nodes tmp
    
  end
  
  if nargin<4
    failure_percent = 0;
  end
  if nargin<3
    n_beacons = 0;
  end
  if n_beacons == 0
    n_beacons = ceil(sqrt(n_nodes));
  end
  k = min(n_beacons, 10); %routing beacon for bvr
  

  if length(topology_fn)==0
    width = sqrt(n_nodes/density);
    [D coord] = simulate_topology(n_nodes, width, width);
  end
  
  tmp = randperm(n_nodes);
  beacons = tmp(1:n_beacons);
  
  %tmp = randperm(n_nodes); %uncomment this line if failure in beacons are allowed
  failure_nodes = tmp(n_nodes-floor(n_nodes*failure_percent)+1:n_nodes);
  if length(failure_nodes)<=20
    failure_nodes
  end
  
  good_nodes = setdiff(1:n_nodes, failure_nodes);
  random_pairs_idx = ceil(rand(n_random_pairs, 2)*length(good_nodes));
  random_pairs = reshape(good_nodes(random_pairs_idx(:)), n_random_pairs, ...
			 2);
  %%clear random_pairs_idx good_nodes tmp
  save topo-tmp.mat D coord density k n_nodes n_beacons failure_nodes ...
      random_pairs n_random_pairs failure_percent beacons
end
fprintf(['# of nodes:%d, density: %.2f, # of test pairs: %d\n# of' ...
	 ' beacons: %d, # of routing beacons: %d\n'], n_nodes, density, ...
	n_random_pairs, n_beacons, k);
t = sprintf( '%d nodes with density %.2f', n_nodes, ...
	     density);
if n_nodes<=200 & show_all_figures
  figure(1);
  visualize(D, coord, beacons,failure_nodes, n_nodes<=50, 1);
  title(t);
else
  fprintf('n_nodes>200 or forbidden to show figures, omit visualization\n');
end

if run_bvr
  %[routing_path, routing_dist_bvr, load_bvr, storage_bvr, dht_bvr, failure_bvr]=bvr(D, n_beacons, ...
	%k , beacons, random_pairs, 1, failure_nodes, 0, 0);
  [routing_path, routing_dist_bvr_sf, load_bvr_sf, storage_bvr, dht_bvr, failure_bvr_sf]=bvr(D, n_beacons, ...
	k , beacons, random_pairs, 1, failure_nodes, 0, 1);
  failure_bvr_sf_real = length(find(routing_dist_bvr_sf==Inf))/ ...
	  n_random_pairs;
else
  load_bvr = 0;
  dht_bvr = -1; 
end
  
if run_bvr_edist  
  [routing_path, routing_dist_bvr2, load_bvr2, storage_bvr2, dht_bvr2, failure_bvr2]=bvr(D, n_beacons, ...
	k , beacons, random_pairs, 0, failure_nodes);
  failure_bvr2_real = length(find(routing_dist_bvr2==Inf))/n_random_pairs;
else
  load_bvr2=0;
  dht_bvr2 = -1;
  failure_bvr2 = 1;
  failure_bvr2_real = 1;
  storage_bvr2 = 0;
  dht_bvr2 = 0;
  
end


[routing_dist_cr, load_cr,storage_cr, dht_cr, cluster_cr, dv_packets_cr, dv_entries_cr] = ...
  compactrouting(D, beacons, random_pairs,failure_nodes, 2, coord);%distance vector
failure_cr = length(find(routing_dist_cr==Inf))/n_random_pairs;

if 0
  [routing_dist_cr, load_cr,storage_cr, dht_cr, cluster_cr] = ...
      compactrouting(D, beacons, random_pairs,failure_nodes, 0);%default, proactive one
  failure_cr = length(find(routing_dist_cr==Inf))/n_random_pairs;
end
if run_cr_reactive
  [routing_dist_cr2, load_cr2,storage_cr2, dht_cr2, cluster_cr2, dv_packets_cr2, dv_entries_cr2] = ...
      compactrouting(D, beacons, random_pairs,failure_nodes, 1, coord); %reactive
  failure_cr2 = length(find(routing_dist_cr2==Inf))/n_random_pairs;
else
  routing_dist_cr2 = routing_dist_cr;
  load_cr2 = load_cr;
  storage_cr2 = storage_cr;
  dht_cr2 = dht_cr;
  load_cr2 = 0;
  failure_cr2 = -1;
  cluster_cr2 = [-1];
end
%fprintf('calculating routing distance..\n');
%[routing_dist_cr, pairhop] = routing_distance(D, nexthop_cr);
%[routing_dist_cr, pairhop] = routing_distance(D, nexthop_cr, random_pairs);

if run_optimal
  neighbor_with_failure =D;
  for i=failure_nodes
    neighbor_with_failure(:,i)=zeros(n_nodes,1);
    neighbor_with_failure(i,:)=zeros(1,n_nodes);
  end
  [optimal_dist optimal_nexthop] = allpair_dijkstra(neighbor_with_failure, ...
						    1:n_nodes, 1, 'optimal-cache.mat');
  pair_indices = random_pairs(:,1)+(random_pairs(:,2)-1)*n_nodes;
  optimal_dist = optimal_dist(pair_indices);
  failure_opt = length(find(optimal_dist==Inf))/n_random_pairs;
else
  optimal_dist = [];
  optimal_nexthop = [];
end


%[routing_dist_bvr routing_dist_cr optimal_dist]
if show_all_figures
  figure(2);
  hold on
end
if run_optimal
  stretch_cr = plot_routing_stretch(routing_dist_cr,optimal_dist,'g',  stretch_percentile, show_all_figures);
  plot_routing_stretch(routing_dist_cr2,optimal_dist,'r.', stretch_percentile, show_all_figures);
  if run_bvr
    stretch_bvr = plot_routing_stretch(routing_dist_bvr_sf,optimal_dist,'b', stretch_percentile, ...
						 show_all_figures);
  end
  if run_bvr_edist
    plot_routing_stretch(routing_dist_bvr2,optimal_dist,'b.', stretch_percentile, show_all_figures);
  end
end
if show_all_figures
  xlim([1 4.1])
  xlabel('routing stretch');
  ylabel('CDF');
  if run_bvr
    legend('cr','cr reactive','bvr','bvr edist',0);
  else
    legend('cr','cr reactive', 0);
  end
  
  title(t);
end

if show_all_figures
  figure(3);
  hold on
  plot_abs_error(routing_dist_cr, optimal_dist,'g');
  plot_abs_error(routing_dist_cr2, optimal_dist,'r.');

  if run_bvr
    plot_abs_error(routing_dist_bvr_sf, optimal_dist,'b');
  end
  if run_bvr_edist
    plot_abs_error(routing_dist_bvr2, optimal_dist,'b.');
  end
  ylabel('CDF');
  xlabel('abs routing error');
  if run_bvr
    legend('cr','cr reactive','bvr','bvr edist',0);
  else
    legend('cr','cr reactive', 0);
  end
  
  title(t);
  fprintf('Done.\n');
end

%%%compare routing load
if run_optimal
  load_optimal = get_nexthop_load(optimal_nexthop, random_pairs);
else
  load_optimal = zeros(length(random_pairs),1);
end

%load_cr = get_nexthop_load(nexthop_cr, random_pairs);
fprintf('==========\n');
fprintf('dht packets: cr:%d \tcr-reactive %d \tbvr %d \tbvr-edist %d\n\n', dht_cr,dht_cr2, ...
	dht_bvr,dht_bvr2);

fprintf('CR mean/max cluster size: %.1f/%d, %.1f/%d\n\n', ...
	mean(cluster_cr), max(cluster_cr), mean(cluster_cr2), ...
	max(cluster_cr2));

fprintf('\nCR failure rate: %.2f%%, CR2 %.2f%%\n', ...
	failure_cr*100, failure_cr2*100);
if run_bvr
  fprintf('BVR failure rate: sf:%.2f%%\n', failure_bvr_sf*100);
  fprintf('BVR sf real failure rate: %.2f%%, %.2f%%\n', failure_bvr_sf_real*100, failure_bvr2_real*100);
end
if run_optimal
  fprintf('Optimal failure rate: %.2f%%\n',failure_opt*100);
end
mmmm(storage_cr,'storage cr',1);
mmmm(storage_cr2,'storage cr-2',0);
if run_bvr
  mmmm(storage_bvr,'storage bvr',0);
  mmmm(storage_bvr2,'storage bvr-edist',0);
end
fprintf('\n');


mmmm(load_cr,'load cr',0);
mmmm(load_cr2,'load cr-2',0);
if run_bvr
  mmmm(load_bvr_sf,'load bvr',0);
  mmmm(load_bvr2,'load bvr-edist',0);
end
if run_optimal
  mmmm(load_optimal,'load opt',0);
end

%fprintf('max load: optimal, cr, bvr = %d,%d,%d\n',max(load_optimal), ...
%	max(load_cr),max(load_bvr));
%fprintf('mean load: optimal, cr, bvr = %.2f,%.2f,%.2f\n',mean(load_optimal), ...
%	mean(load_cr),mean(load_bvr));
if show_all_figures
  figure(4);
  plot_cumu_matrix(load_optimal,'k');
  hold on
  plot_cumu_matrix(load_cr,'g');
  plot_cumu_matrix(load_cr2,'r.');
  if run_bvr
    plot_cumu_matrix(load_bvr_sf,'b');
  end
  if run_bvr_edist
    plot_cumu_matrix(load_bvr2,'b.');
  end
  legend('optimal','cr','cr reactive','bvr','bvr edist',0);
  title(t);
  ylabel('CDF');
  xlabel('node load');
end

res = struct('routing_dist_cr',routing_dist_cr, ...
	     'load_cr', load_cr, 'beacons',beacons, ...
	     'failure_bvr',failure_bvr_sf, ...
	     'failure_cr', failure_cr, ...
	     'optimal_dist',optimal_dist, ...
	     'cluster_cr', cluster_cr, 'routing_dist_bvr_sf', ...
	     routing_dist_bvr_sf, 'load_bvr_sf', load_bvr_sf, ...
	     'dv_packets_cr', dv_packets_cr, 'dv_entries_cr', dv_entries_cr);


function mmmm(m, str, header)
%print mean, max, min, median
m = m(:);
if header
  fprintf('%20s\tmean \t max \t min \t median\n','');
end
fprintf('%20s\t%.1f \t %.1f \t %.1f \t %.1f\n', str, mean(m),max(m), ...
	min(m),median(m));
return