function show_graph( nodelist, xc_len, threshold, num_std, fig )
% nodelist is a 5-digit integer [row1|row2|row3|row4|row5] that plots rowX
% if rowX is non-zero; or a list of nodes [ node1 node2 ... ] e.g. 1:64
% if not specified, all 22 active nodes [25..64] will be shown.
% can pass [] or omit other params to use default values.

if (~exist('fig') || isempty(fig) )
    fig = 580;
end
figure(fig); clf; set(fig,'DoubleBuffer','on') %Flash-free rendering

if (~exist('num_std') || isempty(num_std) )
    num_std = 5;
end
load( sprintf( 'xc%d_data', num_std ) );

if (~exist('threshold') || isempty(threshold) )
    threshold = 4;
end

if (~exist('xc_len') || isempty(xc_len) )
    xc_len = 50;
end
xc1000_len = xc_len/10;

noderow{1} = [ 25 26 27 30 31 ];
noderow{2} = [ 35 36 ];
noderow{3} = [ 41 43 44 46 ];
noderow{4} = [ 49 50 51 52 54 ];
noderow{5} = [ 57 58 59 60 63 64 ];
if ( exist('nodelist') && length(nodelist) )    % exists and non-empty
    if ( length(nodelist) == 1 )
        nodes = [];
        for i = 1:5
            if ( mod( floor(nodelist/10^(5-i)), 10) )
                nodes = [ nodes noderow{i} ];
            end
        end
    else
        nodes = nodelist;
    end
else
    nodes = [ 25 26 27 30 31 35 36 41 43 44 46 49 50 51 52 54 57 58 59 60 63 64 ];
end

channel_pairs = nchoosek( nodes, 2 );
channel_pairs = sortrows( channel_pairs );

edges = zeros(64);
for pair = 1:size(channel_pairs,1)
        ChA = channel_pairs( pair, 1 );
        ChB = channel_pairs( pair, 2 );

        my_xc = xc1000{ChA,ChB};
        my_xc_len = floor( length(my_xc)/2 );
        my_xc = my_xc( my_xc_len+1-xc1000_len : my_xc_len+1+xc1000_len );
        my_xc( xc1000_len : xc1000_len+2 ) = 0;
        
        [maxL,maxLi] = max( my_xc(1:xc1000_len) );
        [maxR,maxRi] = max( my_xc(xc1000_len+2:end) );
        if ( maxL > threshold )
            edges( ChB, ChA ) = maxL;
        end
        if ( maxR > threshold )
            edges( ChA, ChB ) = maxR;
        end
end


rand( 'state', 0 );
warning off MATLAB:divideByZero;
LTable = 1./eye(64) - 1;            % put zeros on diagonal and inf elsewhere
i = 0;
excised_nodes = [];
for node = nodes
    i = i + 1;
    if ( isempty( find( [ edges(node,:) edges(:,node)' ] ) ) )  % doesn't connect to anything, nothing connects to it
        nodes = [ nodes(1:i-1) nodes(i+1:end) ];    % excise node out of nodes
        i = i - 1;
        excised_nodes(end+1) = node;
        continue;                                   % don't plot
    end
    y = 8 - floor((node-1)/8) + .5*(rand-.5);
    x = mod(node-1,8) + 1 + .5*(rand-.5);
    node_coordinates(node,:) = [x y];
    plot( x, y, 'ko', 'linewidth', 2, 'markersize', 20 ); hold on;
    tx = text( x, y, sprintf('%d', node) );
    set( tx, 'fontsize', 10, 'fontweight', 'bold', 'horizontalalignment', 'center' );
    LTable(node,:) = ShortestPath( node, 1, LTable(node,:), edges );
end

running_length=0; running_neighbors=0; running_cluster=0;
i = 0;
for node = nodes
    i = i + 1;
    neighbors = find( LTable(node,:) < inf );   % includes acquaintances (non-immediate neighbors)
    neighbors = neighbors( find( LTable(node,neighbors) > 0 ) );    % exclude the node itself (with length 0)
    num_neighbors = length( neighbors );
    total_length = sum( LTable(node,neighbors) );
    avg_length = total_length / num_neighbors;
    running_length = running_length + total_length;
    running_neighbors = running_neighbors + num_neighbors;
    disp_string = strcat( 'Node %2d / [%2d] neighbors:', repmat( ' %2d', 1, num_neighbors ) );
    disp( sprintf( disp_string, node, num_neighbors, neighbors ) );
    
    degree(i) = length( find( edges(node,:) ) );        
    disp_string = strcat( '        / path length   :', repmat( ' %2d', 1, num_neighbors ), '   <%.2f> / degree %d' );
    disp( sprintf( disp_string, LTable(node,neighbors), avg_length, degree(i) ) );

    neighbors = find( edges(node,:) );          % only adjacent neighbors
    num_neighbors = length( neighbors );
    if ( num_neighbors > 1 )
        pairs = nchoosek( neighbors, 2 );
        denom = num_neighbors * (num_neighbors-1);
    else
        pairs = [];
        denom = 1;
    end
    
    total_edges = 0;
    for pair = 1:size(pairs,1)
        ChA = pairs( pair, 1 );
        ChB = pairs( pair, 2 );
        if ( edges( ChA, ChB ) )
            total_edges = total_edges + 1;
        end
        if ( edges( ChB, ChA ) )
            total_edges = total_edges + 1;
        end        
    end
    running_cluster = running_cluster + total_edges/denom;

    disp( sprintf( '        / C(%d) = %d / %d = %.2f', node, total_edges, denom, total_edges/denom ) );
end


disp( sprintf( '\nTotal # nodes = %d', length(nodes) ) );
disp_string = strcat( 'Excised nodes = ', repmat( ' %2d', 1, length(excised_nodes) ) );     % no neighbors AND not any node's neighbor
disp( sprintf( disp_string, excised_nodes ) );
L = running_length / running_neighbors;
disp( sprintf( '<length>, L   = %d / %d  = %.3f', running_length, running_neighbors, L ) );
C = running_cluster / length(nodes);
disp( sprintf( '<C(i)>, C     = %.2f / %d = %.3f       1/N = %.3f      #times larger than 1/N = %.2f', running_cluster, length(nodes), C, ...
    1/length(nodes), C*length(nodes) ) );

axis off;
tt = title( sprintf( 'Edge Connectivity Graph:  xc\\_len=%dms  xc\\_threshold=%.1f  num\\_std=%.1f', xc_len, threshold, num_std ) );
set (tt, 'fontsize', 12, 'fontweight', 'bold' );
AX = axis;
tx = text( 0.1*AX(1)+0.9*AX(2), 0.4*AX(3)+0.6*AX(4), ...
           sprintf( 'L = %.3f\nC = %.3f\n# nodes, N = %d\nC > 1/N by %.1f', L, C, length(nodes), C*length(nodes) ) );
set (tx, 'fontsize', 12, 'fontweight', 'bold' );
drawnow;

save( sprintf( 'graph%d_data', num_std ), 'edges', 'LTable', 'nodes', 'node_coordinates', 'L', 'C' );
%disp( 'Graph saved. Hit a key to continue: ' );
%pause;

colormap('default');
default_colormap = colormap;
r = 0.2;
[ChA,ChB,val] = find( edges );
max_val = max(val);
min_val = min(val);
for edgenum = 1:length(ChA)
    u = node_coordinates( ChA(edgenum), : );
    v = node_coordinates( ChB(edgenum), : );
    w = v - u;
    w = w / sqrt( w(1)^2 + w(2)^2 );
    u = u + r*w;
    v = v - r*w;
    color = round( (val(edgenum)-min_val)/(max_val-min_val) * 63 + 0.5);
    h = arrowline( [u(1) v(1)], [u(2) v(2)], 'color', default_colormap(color,:), 'arrowsize', 200 );
%    disp( sprintf( '%2d->%2d : %.2f', ChA(edgenum), ChB(edgenum), val(edgenum) ) );    
    drawnow;
end


figure(fig+1); clf;
hist( degree ); hold on;

[myY,myX] = hist( degree );
Y=[]; X=[];
for i = find( myY )
    X(end+1) = myX(i);
    Y(end+1) = myY(i);
end

[P,S] = polyfit( log10(X), log10(Y), 1 );
fitY = 10^P(2) * myX.^P(1);
plot( myX, fitY, 'go-', 'linewidth', 2', 'markersize', 6 );

xlabel( 'Degree, k', 'fontsize', 12, 'fontweight', 'bold' );
ylabel( 'Number of occurences', 'fontsize', 12, 'fontweight', 'bold' );
tt = title( sprintf( 'Histogram of Degree Distribution  /  xc\\_len=%dms  xc\\_threshold=%.1f  num\\_std=%.1f', xc_len, threshold, num_std ) );
set (tt, 'fontsize', 12, 'fontweight', 'bold' );
AX = axis;
tx = text( 0.25*AX(1)+0.75*AX(2), 0.05*AX(3)+0.95*AX(4), ...
           sprintf( '# nodes = %d\n\\gamma = %.2f', length(nodes), P(1) ) );
set (tx, 'fontsize', 12, 'fontweight', 'bold' );       

return;