function B = btermatch(G, graphname)
% BTERMATCH : Match a BTER graph to a given graph G
%
% B = btermatch(G);
% B = btermatch(G, graphname);
%
% Script from Feastpack web pages on matching BTER to a given graph G

if nargin < 2
  graphname = 'input graph';
end;

% Some preprocessing may be needed to make the graph simple.

G = G + G'; %<- Make the graph undirected
G = spones(G); %<- Make the graph unweighted
% G = spdiags(zeros(size(G,1),1),0,G); %<- Remove self-loops
G = G - diag(diag(G));


% We use an example graph from the SNAP collection which has already been preprocessed and is included in the distribution.
% 
% load('data/ca-AstroPh-simple');

nnodes = size(G,1);
nedges = nnz(G)/2;
fprintf('\nCreating BTER graph to match input graph\n\n');
fprintf('Graph name: %s\n', graphname);
fprintf('Number of nodes: %d\n', nnodes);
fprintf('Number of edges: %d\n', nedges);

% Graph name: ca-AstroPh-simple
% Number of nodes: 18772
% Number of edges: 198050

% Step 2: Compute the degree distribution

% Calculate the number of nodes per degree.

nd = accumarray(nonzeros(sum(G,2)),1);
maxdegree = find(nd>0,1,'last');
fprintf('Maximum degree: %d\n', maxdegree);

% Maximum degree: 504

% Step 3: Calculate the clustering coefficient per degree

% This can be extremely expensive, especially for large graphs. We
% provide C code and accompanying mex interfaces to do both exact and
% estimated calculations. The matlab function ccperdeg is an interface
% to both; see below for details. The ccperdeg function returns the
% clustering coefficient per degree (ccd) and the global clustering
% coefficient (gcc). In the second case, the value of 'nsamples' is
% the number of sample wedges per bin.

[ccd,gcc] = ccperdeg(G);
fprintf('Global clustering coefficient: %.2f\n', gcc);

% Global clustering coefficient: 0.32

% Step 4: Generate BTER model

% BTER creates two edges lists, corresponding to Phase 1 and Phase 2 edges, respectively.

fprintf('\nRunning BTER...\n');
t1=tic;
[E1,E2] = bter(nd,ccd);
toc(t1)
fprintf('Number of edges created by BTER: %d\n', size(E1,1) + size(E2,1));

% Running BTER...
% --- BTER HPC Set-up ---
% Desired number of nodes: 18771
% Desired number of edges: 198050
% Multiplier to degree-1 nodes: 1
% Maximum degree: 504
% Number of groups: 76
% Number of blocks: 2095
% Phase 1 total weight: 303303
% Phase 2 total weight: 52313
% Time for setup (sec): 0.02
% Determined phase for 355616 edges in 0.018912 seconds
% Sampled 303534 edges for phase 1 in 0.185769 seconds
% Sampled 52082 edges for phase 2 in 0.058720 seconds
% Removed 35 loops in 0.001642 seconds
% --- BTER HPC Complete ---
% Elapsed time is 0.319444 seconds.
% Number of edges created by BTER: 355581

% Step 5: Create sparse matrix representation of BTER graph

% Converting the edge list to an adjacency matrix is often the most
% expensive part in MATLAB. We separate this step because the edge
% list may be sufficient for many applications.

fprintf('Turning edge list into adjacency matrix (including dedup)...\n');
t2=tic;
G_bter = bter_edges2graph(E1,E2);
toc(t2);
B = G_bter;
fprintf('Number of edges in dedup''d graph: %d\n', nnz(G)/2);

% Turning edge list into adjacency matrix (including dedup)...
% Elapsed time is 0.770594 seconds.
% Number of edges in dedup'd graph: 198050

% Step 6: Calculate metrics on BTER graph and compare to original

% Number of nodes and edges
nnodes_bter = size(G_bter,1);
nedges_bter = nnz(G_bter)/2;
fprintf('\nComparative Results: %s versus BTER\n', graphname);
fprintf('Number of nodes: %d versus %d\n', nnodes, nnodes_bter);
fprintf('Number of edges: %d versus %d\n', nedges, nedges_bter);

% Degree disttribution
nd_bter = accumarray(nonzeros(sum(G_bter,2)),1);
maxdegree_bter = find(nd_bter>0,1,'last');
fprintf('Maximum degree: %d versus %d\n', maxdegree, maxdegree_bter);

% Clustering coefficients
[ccd_bter,gcc_bter] = ccperdeg(G_bter);
fprintf('Global clustering coefficient: %.2f versus %.2f\n\n', gcc, gcc_bter);

% Comparative Results: ca-AstroPh-simple versus BTER
% Number of nodes: 18772 versus 18771
% Number of edges: 198050 versus 195795
% Maximum degree: 504 versus 423
% Global clustering coefficient: 0.32 versus 0.33

% Step 7: Visualize the data

% Degree distribution
figure;
loglog(1:length(nd),nd,'ro',1:length(nd_bter),nd_bter,'b*');
title('Degree distribution');
xlabel('Degree');
ylabel('Frequency');
legend(graphname, 'BTER');

% Clustering coefficients
figure;
semilogx(1:length(ccd),ccd,'ro',1:length(ccd_bter),ccd_bter,'b*');
title('Clutering Coefficients');
xlabel('Degree');
ylabel('Mean CC');
legend(graphname, 'BTER');

% Step 7 (alternate): Visualize the results with binning

% The data can be extremely crowded, especially for the higher degrees.
% We also include the bindata function that can be used to logrithmically
% bin the data. Type 'help bindata' for full instructions. We make
% the following parameter choices here:
% 
%     tau = 10 : The first 10 bins are singleton bins.
%     omega = 1.5 : Subsequent bins grow in size by a factor of 1.5.
%     ybinfun = @mean : Take the average of all entries in the bin.
%     prebin = true : Ensures every entry in the bin has an assigned value, 
%     even if its zero. Especially significant for last bin which may not 
%     be completely specified by the inputs.

binparams = {'tau', 10,'omega', 1.5, 'ybinfun', @mean, 'prebin', false};

% In the binned version, the degree distibution can be non-integer
% and less than one because its the average value for the bin.

% Degree distribution
figure;
[xx1,yy1] = bindata([], nd, binparams{:});
[xx2,yy2] = bindata([], nd_bter, binparams{:});
loglog(xx1,yy1,'ro',xx2,yy2,'b*');
title('Degree distribution (logrithmically binned with means)');
xlabel('Degree');
ylabel('Frequency');
legend(graphname, 'BTER');

% Clustering coefficients

figure;
[xx1,yy1] = bindata([], ccd, binparams{:});
[xx2,yy2] = bindata([], ccd_bter, binparams{:});
semilogx(xx1,yy1,'ro',xx2,yy2,'b*');
title('Clutering Coefficients (logrithmically binned with means)');
xlabel('Degree');
ylabel('Mean CC');
legend(graphname, 'BTER');


% "Blowup" option in BTER for better match to degree-1 vertices

% In many cases, BTER underestimates the number of degree-1 nodes.
% This is a well-understood problem, and we use the correction
% prescribed by Durak et al. (arXiv:1210.5288, 2012). We can see that
% now the number of degree-1 vertices is much more accurate.

[E1_alt,E2_alt] = bter(nd,ccd,'verbose',false,'blowup',10);
G_alt = bter_edges2graph(E1_alt, E2_alt);
nd_alt = accumarray(nonzeros(sum(G_alt,2)),1);
figure;
loglog(1:length(nd),nd,'ro',1:length(nd_bter),nd_bter,'b*');
hold on;
loglog(1:length(nd_alt),nd_alt,'g^','Color', [0 .7 .2]);
title('Degree distribution');
xlabel('Degree');
ylabel('Frequency');
legend(graphname, 'BTER', 'BTER with Blowup');
