%-------------------------------------
% Lecture 12. 
% November 26, 2010.
% Complexity in networs and brain


% Hilbert curve simulation
x=0; y=0;
newx=.5*[-.5+y -.5+x .5+x  .5-y]
newy=.5*[-.5+x  .5+y .5+y -.5-x]

figure;plot(newx,newy,'-')

x=newx;
y=newy;

newx=.5*[-.5+y -.5+x .5+x  .5-y]
newy=.5*[-.5+x  .5+y .5+y -.5-x]

figure;plot(newx,newy,'-')


%function [x,y] = hilbert(n)
%if n<=0
%    x=0;
%    y=0;
%else
%    [xo,yo]=hilbert(n-1); 
%    x=.5*[-.5+yo -.5+xo .5+xo  .5-yo];
%    y=.5*[-.5+xo  .5+yo .5+yo -.5-xo];
%end

n=5;
[x,y]=hilbert(n);
figure;line(x,y)

%--------------------------------------
% Fractal dimension (FD) of point cloud data

% 1D pattern with noise
x=1:100;
y=x + normrnd(0,5,1,100);
figure; plot(x,y,'.');


% 2D pattern with noise
a=1:100;
x=a + normrnd(0,20,1,100);
y=a + normrnd(0,20,1,100);
figure; plot(x,y,'.');

% compute the pairwise distance

D=zeros(100,100);
for i=1:100
    D(i,:)=sqrt((x-kron(ones(1,100),x(i))).^2 + (y-kron(ones(1,100),y(i))).^2);
end;

figure; imagesc(D); colorbar


% correlation dimension
maxk=floor(max(max(D)));
D=reshape(D,100^2,1);
count=[];
for k=1:maxk
    count(k)= length(find(D<k))/2;
end;
count=count/maxk^2;

figure; plot(log(1:maxk),log(count), '.')
b = regress(log(count)', [ones(maxk,1) log(1:maxk)'])

%--------------------------
% box counting method
C=imread('fractal.gif');
[N,R]=boxcount(C,'plot');

%---------------------------------------
% box counting method applied to network

% fiber bundle visualization
load SL.mat
for i=1:30:10000 % subsampling 1 for every 30 tracts
    tract=SL{i}'; 
    if size(tract,2)> 10
        hold on
        plot3(tract(1,:),tract(2,:),tract(3,:),'b')
    end;
end;

% obtain end points of tracts
endp=get_endp_tracts(SL,10);

hold on; plot3(endp(:,1,1),endp(:,1,2), endp(:,1,3),'.r')
hold on; plot3(endp(:,2,1),endp(:,2,2), endp(:,2,3),'.r')
set(gcf,'Color','white','InvertHardcopy','off');

% epsilon-neighbor graph construction
[adj, prob, vertices]=tract2graph(SL, 6);
figure_dtigraph(adj,vertices);
figure; imagesc(adj)
figure; [N,R]=boxcount(adj,'plot');

k=length(N);
b = regress(log(N)', [ones(k,1) log(R)'])
% FD=1.3385

% without symmetry
n=size(adj,1);
T=zeros(n,n);
for i=1:n
  T(i,i:end)=adj(i,i:end);
end;
figure; [N,R]=boxcount(T,'plot');

k=length(N);
b = regress(log(N)', [ones(k,1) log(R)'])
% FD=1.2341

%-----------------------
% population study on FD
% this is only for class demonstration
% data is not uploaded into webspace

% id=textread('REPORT.tracts.oct.2009.txt','%s')
% 
% figure;
% FD=[];
% for i=1:14 %control subjects
%     i
%     SL=get_streamlines(id{i},[1.5 1.75 2.25]);
%     [adj, prob, vertices]=tract2graph(SL,10);
%     [N,R]=boxcount(adj,'slope');
%     %[N,R]=boxcount(adj,'plot');
%     hold on;
%     k=length(N);
%     b = regress(log(N)', [ones(k,1) log(R)']);
%     FD(i)=-1*b(2);
% end;
% %FD = 1.3434 +- 0.0513
% set(gcf,'Color','w');
% print -r300 -dtiff FD_control_slope

%---------------------------
% random graphs

figure;
FD=[];
for j=1:14
    % constuct random adjacency matrix of size 200 x 200
    adj=binornd(1, 0.99, 200,200);
    n=size(adj,1);
    T=zeros(n,n);
    for i=1:n
        T(i,i:end)=adj(i,i:end);
    end;
    T=T+T';
    T=T-diag(diag(T));
    [N,R]=boxcount(T,'plot');
    hold on;
end;

set(gcf,'Color','w');
%print -r300 -dtiff FD_binary0.99_slope


