function back = dd( A , type , quiet)
%DD  determines and displays the degree distribution of a network
%  DD(A) plots the degree distribution of a network
%
%  DD(A,'type')  where 'type' may be 'er' for Erdos-Renyi, 'sw'
%     for Small World or 'sf' for Scale Free specialised recursion
%
%  Last change:  10/2/2005 - Florian Knorn

set(0,'defaulttextinterpreter','none'); warning off;
figure(1); hold on;
lw = 1; % linewidth for plots im dokument: 1 (!)
pw = 1; % width of the + markers
plotfit = 1;

if nargin < 1
	error('Please input adjacency matrix !');
else
	if nargin == 1
		type = 'none';
	end
	n = size(A,1);
	if n ~= size(A,2)
		error('Please input SQUARE matrix !');
	end
end

if nargin ~= 3
	quiet = false;
end

degs = sum(full(A),2);	% orinial degrees of each node

[d,i] = max(degs);		% d - max. degree, i - the node index with this max degree
if ~quiet
	disp(sprintf('Highest degree: Node %i with degree %i',i,d));
end

if strcmp(type,'sf')
	degree = 1:d;
else
	degree = 0:d;
end
count = histc(degs,degree);% sort degrees into bins 0 ... d
count = count'/n;			% count -> fractions

% draw data with values in % in a linlin plot
if ~strcmp(type,'sf')% subplot(2,1,2); end
	plot(degree,100*count,'k+','linewidth',pw); hold on;
end

% just some some fancy labeling
if max(degree) < 60 && max(degree) > 20 && ~strcmp(type,'sf')
	labelz = char('0');
	for i = 1:2:d
		if mod(i,2)
			labelz = char(labelz , int2str(i));
		end
	end
	if mod(d,2) % uneven number
		set(gca,'XLim',[0 d+1],'XTick',[0 1:2:d],'XTickLabel',labelz);
	else % even number -> add last label
		set(gca,'XLim',[0 d+1],'XTick',[0 1:2:d d],'XTickLabel',char(labelz,int2str(d)));
	end
else
	set(gca,'XLim',[0 max(degree)]);
end

xlabel('Node Degree \incarr'); ylabel('Fraction of nodes in \% \incarr');


if strcmp(type,'er') && plotfit % ================= ERDÖS RENYI =====================================

	% first: remove all entries containing zeros
	combined = [ degree ; count ];

	[i,j] = find(combined(2,:) == 0);	% find the zeros
	combined(:,j) = [];				% and kick them out

	nzdegree = combined(1,:); nzcount = combined(2,:);
	% second: determine and plot parameters by nonlinear fitting
	if n < 5100
		% use binomial distribution
		global Nm1 binom; binom = 1; Nm1 = n-1;
		param0 = avk(A)/n;		% parameter guess for p, to get fitting started
		param = nlinfit(nzdegree,nzcount,@er_reg_model,param0);	% nonlinear fitting

		if ~quiet		
			disp(sprintf('Fitted binomial distribution resulted in a guessed p = %f',param));
		end

		% for plotting the fitted model:
		gk = 0:degree(end); % x-coordinates

		gcount = er_reg_model(param,gk); % y-coord
	else
		% use poisson distribution
		global binom; binom = 0;
		param0 = avk(A);		% parameter guess for <k>, to get fitting started
		param = nlinfit(nzdegree,nzcount,@er_reg_model,param0);	% nonlinear fitting
		if ~quiet
			disp(sprintf('Fitted Poisson distribution resulted in a guessed <k> = %f',param));
		end

		% for plotting the fitted model:
		gk = 0:degree(end); % x-coordinates
		gcount = er_reg_model(param,gk);	 % y-coord
	end

	% plot fitted curve
	plot(gk,gcount*100,'r--','linewidth',lw); hold on;

elseif strcmp(type,'sw') && plotfit % ================= SMALL WORLD =================================

	% remove all entries containing zeros
	combined = [ degree(2:end) ; count(2:end) ]; % combine vectors
	[i,j] = find(combined == 0);	% find the zeros-count columns...
	combined(:,j) = [];				% ... and remove them

	% little things...
	nzdegree = combined(1,:); nzcount = combined(2,:); global K; leg = 0;

	% find a guess for K by picking the two closest pair integers to avk(A)
	K = avk(A);

	% guess parameters, to get nonlinear fitting started param(1) = p, param(2) = scaling
	param0 = [0.5 n];

	% now start the fitting
	[param,R1,J1] = nlinfit(nzdegree,nzcount,@sw_reg_model,param0);

	if ~quiet
		disp(sprintf('Fitting using K = %i succeeded, resulting in p = %.2f',K,param(1)));
	end

	gk = min(nzdegree):d; gcount = sw_reg_model(param,gk);
	plot(gk,gcount,'r--','linewidth',lw); % in loglog

	% subplot(2,1,2);
	param = [0.01 3000]; % cheating ;-)
	gk = 0:d; gcount = sw_reg_model(param,gk);
	plot(gk,gcount*100,'r--','linewidth',lw); % in normal plot below

elseif strcmp(type,'sf') % ================= SCALE FREE ==================================
	han1 = figure(1); ch = get(gcf,'Children');
	if length(ch) == 2
		han2 = ch(2);
		set(han1,'CurrentAxes',ch(1));
	else
		han2 = get(han1,'CurrentAxes');
		pos = get(han2,'Pos');
		x1 = pos(1); w1 = pos(3);	y1 = pos(2); h1 = pos(4);
		w2 = 0.45; 	h2 = 0.45; % width and height of the inset
		x2 = x1 + (w1-w2); y2 = y1 + (h1-h2);
		han3 = axes('Pos',[x2 y2 w2 h2]);
	end

	hold on;
	% draw log-log-plot
	plot(degree,count*100,'k+','linewidth',pw); 
	set(gca,'XScale','log','YScale','log','XMinorTicks','on','YMinorTicks','on'); lims = axis;
	
	% first: remove all entries containing zeros
%-----------
	combined = [ degree(1:end) ; count(1:end) ];

	[i,j] = find(combined == 0);	% find the zeros
	combined(:,j) = [];				% and kick them out

	nzdegree = combined(1,:); nzcount = combined(2,:);

	param0 = [max(count) 2 1];		% parameter guess, to get recursion started

	% second: determine parameters by recursion
	param = nlinfit(nzdegree,nzcount,@sf_reg_model,param0);
	back = param(1);

	if plotfit
		% plot fitted curve in loglog
		disp(sprintf('Fitted powerlaw distribution: p(k) = %f (k+%f)^(-%f)',param(1),param(3),param(2)));

		gk = [1:d+10]; gcount = sf_reg_model(param,gk);
	
		loglog(gk,gcount*100,'r--','linewidth',lw);
		set(gca,'XMinorTicks','on','YMinorTicks','on');		
	end % if plotfit
	
	% plot in linlin (big axes)
	set(han1,'CurrentAxes',han2); hold on;
	plot(degree,100*count,'k+','linewidth',pw);
	if plotfit % plot fitted
		gk = [0.5 1:d]; gcount = sf_reg_model(param,gk);
		plot(gk,gcount*100,'r--','linewidth',lw);
	end

end % switch mode

set(gca,'YLim',[0 max(count)*110]);
set(gca,'XMinorTicks','on','YMinorTicks','on');

warning on;