%% Functional connectivity models

%% 1. Sequential A->B->C

% Unit A is input time series
% Unit properties:
%   N_i         Noise imposed on signal
%
% Connection properties:
%   c_ij        (Synaptic) connection strength
%   d_ij        Axonal delay

%% 1.1 Zero noise, zero delay
F_s = 1000;
T = 1:1/F_s:10;    % ten second time series at 1000Hz
T=T';
L = length(T);
F = 30;            % 30 Hz frequency
F_rps = 2*pi*F;    % In radians/sec
d_ab = 0; d_bc = 0; % Delays
td_ab = round(d_ab * F_s);  % Delay in time points
td_bc = round(d_bc * F_s);
c_ab = 1; c_bc = 1; % Synaptic strengths
N_a = 0; %zeros(L,1);
N_b = N_a; N_c = N_a;
Noise_a = N_a * zscore(randn(L,1)) / 3;
Noise_b = N_b * zscore(randn(L,1)) / 3;
Noise_c = N_c * zscore(randn(L,1)) / 3;

Y_a = sin(F_rps * T) + Noise_a;     % Signal for A is sinusoid plus noise

Y_b = [zeros(td_ab,1) c_ab * Y_a(td_ab+1 : end)] + Noise_b;

Y_c = [zeros(td_bc,1) c_bc * Y_b(td_bc+1 : end)] + Noise_c;

%% 1.2. Non-zero, equal noise, zero delay
N_b = 2;
N_c = 2;

Noise_a = N_a * zscore(randn(L,1)) / 3;
Noise_b = N_b * zscore(randn(L,1)) / 3;
Noise_c = N_c * zscore(randn(L,1)) / 3;

Y_a = sin(F_rps * T) + Noise_a;     % Signal for A is sinusoid plus noise

Y_b = [zeros(td_ab,1) c_ab * Y_a(td_ab+1 : end)] + Noise_b;

Y_c = [zeros(td_bc,1) c_bc * Y_b(td_bc+1 : end)] + Noise_c;

%% 2. Vary across noise levels

N_itr = 10;

Ns = [0 : 0.1 : 10];
relB = [0 : 0.1 : 2];

R = zeros(length(Ns),length(relB),3);

for itr = 1 : N_itr
    for i = 1 : length(Ns)
        N_i = Ns(i); 
        for j = 1 : length(relB)
            rb = relB(j);
            Noise_a = N_i * zscore(randn(L,1)) / 3;
            Noise_b = rb*N_i * zscore(randn(L,1)) / 3;
            Noise_c = N_i * zscore(randn(L,1)) / 3;

            Y_a = zscore(sin(F_rps * T) + Noise_a);

            Y_b = zscore([zeros(td_ab,1) c_ab * Y_a(td_ab+1 : end)] + Noise_b);

            Y_c = zscore([zeros(td_bc,1) c_bc * Y_b(td_bc+1 : end)] + Noise_c);

            R(i,j,1) = R(i,j,1) + corr(Y_a,Y_b);
            R(i,j,2) = R(i,j,2) + corr(Y_b,Y_c);
            R(i,j,3) = R(i,j,3) + corr(Y_a,Y_c);
        end
    end
end

R = R / N_itr;

%% Plot results

conns = {'\rho(AB)','\rho(BC)','\rho(AC)'};
names = {'ab','bc','ac'};

for i = 1 : 3
   h=figure;
   set(h,'Color','w');
   hh=imagesc(squeeze(R(:,:,i)));
   n = round(length(relB) / 10);
   set(gca,'XTickLabels',relB(n:n:end));
   n = round(length(Ns) / 10);
   set(gca,'YTickLabels',Ns(n+1:n:end));
   set(gca,'FontSize',12);
   hh=title(conns{i});
   set(hh,'FontSize',18);
   hh=xlabel('Relative Noise (B)');
   set(hh,'FontSize',16);
   hh=ylabel('Relative Noise (Baseline)');
   set(hh,'FontSize',16);
   caxis([0 1]);
   colorbar;
   
   print(sprintf('seq_img_%s.png',names{i}),'-dpng','-r0');
   %exportfig(h,sprintf('seq_img_%s.png',names{i}));
   %export_fig(sprintf('seq_img_%s.png',names{i}),'-png','-a1');
end

relbs = [0.5 1 1.5];

% Plot lines
for j = 1 : length(relbs)
    i = find(relB==relbs(j));
    h=figure;
    set(h,'Color','w');
    M = squeeze(R(:,i,:));
    hh = plot(Ns,M);
    
    set(gca,'FontSize',12);
    
    num = sprintf('%1.1f', relB(i));
    
    hh=title(['w_b = ' num ' \times baseline']);
    set(hh,'FontSize',18);
    hh=xlabel('Relative Noise (Baseline)');
    set(hh,'FontSize',16);
    hh=ylabel('\rho_{ij}');
    set(hh,'FontSize',16);
    
    hh=legend({'\rho(AB)','\rho(BC)','\rho(AC)'});
    set(hh,'FontSize',14);
    
    print(sprintf('seq_lines_wb_%02d.png',relB(i)*10),'-dpng','-r0');
end



%% Parallel network, vary B noise

d_ab = 0; d_bc = 0; d_bd = 0; % Delays
td_ab = round(d_ab * F_s);  % Delay in time points
td_bc = round(d_bc * F_s);
td_bd = round(d_bd * F_s);
c_ab = 1; c_bc = 1; c_bd = 1; % Synaptic strengths

N_itr = 10;

Ns = [0 : 0.1 : 10];
relB = [0 : 0.1 : 2];

R = zeros(length(Ns),length(relB),6);

for itr = 1 : N_itr
    for i = 1 : length(Ns)
        N_i = Ns(i); 
        for j = 1 : length(relB)
            rb = relB(j);
            Noise_a = N_i * zscore(randn(L,1)) / 3;
            Noise_b = rb*N_i * zscore(randn(L,1)) / 3;
            Noise_c = N_i * zscore(randn(L,1)) / 3;
            Noise_d = N_i * zscore(randn(L,1)) / 3;

            Y_a = zscore(sin(F_rps * T) + Noise_a);

            Y_b = zscore([zeros(td_ab,1) c_ab * Y_a(td_ab+1 : end)] + Noise_b);

            Y_c = zscore([zeros(td_bc,1) c_bc * Y_b(td_bc+1 : end)] + Noise_c);
            
            Y_d = zscore([zeros(td_bd,1) c_bc * Y_b(td_bd+1 : end)] + Noise_d);

            R(i,j,1) = R(i,j,1) + corr(Y_a,Y_b);
            R(i,j,2) = R(i,j,2) + corr(Y_b,Y_c);
            R(i,j,3) = R(i,j,3) + corr(Y_a,Y_c);
            R(i,j,4) = R(i,j,4) + corr(Y_a,Y_d);
            R(i,j,5) = R(i,j,5) + corr(Y_b,Y_d);
            R(i,j,6) = R(i,j,6) + corr(Y_c,Y_d);
        end
    end
end

R = R / N_itr;

%% Plots

conns = {'\rho(AB)','\rho(BC)','\rho(CD)'};
names = {'ab','bc','cd'};
idx = [1 2 6];

for i = 1 : 3
   h=figure;
   set(h,'Color','w');
   hh=imagesc(squeeze(R(:,:,idx(i))));
   n = round(length(relB) / 10);
   set(gca,'XTickLabels',relB(n:n:end));
   n = round(length(Ns) / 10);
   set(gca,'YTickLabels',Ns(n+1:n:end));
   set(gca,'FontSize',12);
   hh=title(conns{i});
   set(hh,'FontSize',18);
   hh=xlabel('Relative Noise (B)');
   set(hh,'FontSize',16);
   hh=ylabel('Relative Noise (Baseline)');
   set(hh,'FontSize',16);
   caxis([0 1]);
   colorbar;
   
   print(sprintf('par_img_%s.png',names{i}),'-dpng','-r0');
end

relbs = [0.5 1 1.5];

% Plot lines
for j = 1 : length(relbs)
    i = find(relB==relbs(j));
    h=figure;
    set(h,'Color','w');
    M = squeeze(R(:,i,idx));
    hh = plot(Ns,M);
    
    set(gca,'FontSize',12);
    
    num = sprintf('%1.1f', relB(i));
    
    hh=title(['w_b = ' num ' \times baseline']);
    set(hh,'FontSize',18);
    hh=xlabel('Relative Noise (Baseline)');
    set(hh,'FontSize',16);
    hh=ylabel('\rho_{ij}');
    set(hh,'FontSize',16);
    
    hh=legend(conns);
    set(hh,'FontSize',14);
    
    print(sprintf('par_lines_wb_%02d.png',relB(i)*10),'-dpng','-r0');
end


%% Parallel network, vary C noise

d_ab = 0; d_bc = 0; d_bd = 0; % Delays
td_ab = round(d_ab * F_s);  % Delay in time points
td_bc = round(d_bc * F_s);
td_bd = round(d_bd * F_s);
c_ab = 1; c_bc = 1; c_bd = 1; % Synaptic strengths

N_itr = 10;

Ns = [0 : 0.1 : 10];
relC = [0 : 0.1 : 2];

R = zeros(length(Ns),length(relC),6);

for itr = 1 : N_itr
    for i = 1 : length(Ns)
        N_i = Ns(i); 
        for j = 1 : length(relC)
            rc = relC(j);
            Noise_a = N_i * zscore(randn(L,1)) / 3;
            Noise_b = 1.5*N_i * zscore(randn(L,1)) / 3;
            Noise_c = rc*N_i * zscore(randn(L,1)) / 3;
            Noise_d = rc*N_i * zscore(randn(L,1)) / 3;

            Y_a = zscore(sin(F_rps * T) + Noise_a);

            Y_b = zscore([zeros(td_ab,1) c_ab * Y_a(td_ab+1 : end)] + Noise_b);

            Y_c = zscore([zeros(td_bc,1) c_bc * Y_b(td_bc+1 : end)] + Noise_c);
            
            Y_d = zscore([zeros(td_bd,1) c_bc * Y_b(td_bd+1 : end)] + Noise_d);

            R(i,j,1) = R(i,j,1) + corr(Y_a,Y_b);
            R(i,j,2) = R(i,j,2) + corr(Y_b,Y_c);
            R(i,j,3) = R(i,j,3) + corr(Y_a,Y_c);
            R(i,j,4) = R(i,j,4) + corr(Y_a,Y_d);
            R(i,j,5) = R(i,j,5) + corr(Y_b,Y_d);
            R(i,j,6) = R(i,j,6) + corr(Y_c,Y_d);
        end
    end
end

R = R / N_itr;

%% Plots

conns = {'\rho(AB)','\rho(BC)','\rho(CD)'};
names = {'ab','bc','cd'};
idx = [1 2 6];

for i = 1 : 3
   h=figure;
   set(h,'Color','w');
   hh=imagesc(squeeze(R(:,:,idx(i))));
   n = round(length(relC) / 10);
   set(gca,'XTickLabels',relC(n:n:end));
   n = round(length(Ns) / 10);
   set(gca,'YTickLabels',Ns(n+1:n:end));
   set(gca,'FontSize',12);
   hh=title(conns{i});
   set(hh,'FontSize',18);
   hh=xlabel('Relative Noise (C,D)');
   set(hh,'FontSize',16);
   hh=ylabel('Relative Noise (Baseline)');
   set(hh,'FontSize',16);
   caxis([0 1]);
   colorbar;
   
   print(sprintf('par_img_varc_%s.png',names{i}),'-dpng','-r0');
end

relbs = [0.5 1 1.5];

% Plot lines
for j = 1 : length(relbs)
    i = find(relC==relbs(j));
    h=figure;
    set(h,'Color','w');
    M = squeeze(R(:,i,idx));
    hh = plot(Ns,M);
    ylim([0 1]);
    
    set(gca,'FontSize',12);
    
    num = sprintf('%1.1f', relC(i));
    
    hh=title(['w_{c,d} = ' num ' \times baseline']);
    set(hh,'FontSize',18);
    hh=xlabel('Relative Noise (Baseline)');
    set(hh,'FontSize',16);
    hh=ylabel('\rho_{ij}');
    set(hh,'FontSize',16);
    
    hh=legend(conns);
    set(hh,'FontSize',14);
    
    print(sprintf('par_lines_varc_wb_%02d.png',relC(i)*10),'-dpng','-r0');
end


%% Parallel network, vary C noise, partial correlations

d_ab = 0; d_bc = 0; d_bd = 0; % Delays
td_ab = round(d_ab * F_s);  % Delay in time points
td_bc = round(d_bc * F_s);
td_bd = round(d_bd * F_s);
c_ab = 1; c_bc = 1; c_bd = 1; % Synaptic strengths

N_itr = 10;

Ns = [0 : 0.1 : 10];
relC = [0 : 0.1 : 2];

R = zeros(length(Ns),length(relC),6);

for itr = 1 : N_itr
    for i = 1 : length(Ns)
        N_i = Ns(i); 
        for j = 1 : length(relC)
            rc = relC(j);
            Noise_a = N_i * zscore(randn(L,1)) / 3;
            Noise_b = 1.5*N_i * zscore(randn(L,1)) / 3;
            Noise_c = rc*N_i * zscore(randn(L,1)) / 3;
            Noise_d = rc*N_i * zscore(randn(L,1)) / 3;

            Y_a = zscore(sin(F_rps * T) + Noise_a);

            Y_b = zscore([zeros(td_ab,1) c_ab * Y_a(td_ab+1 : end)] + Noise_b);

            Y_c = zscore([zeros(td_bc,1) c_bc * Y_b(td_bc+1 : end)] + Noise_c);
            
            Y_d = zscore([zeros(td_bd,1) c_bc * Y_b(td_bd+1 : end)] + Noise_d);

            Y = [Y_a Y_b Y_c Y_d];
            
            R_ij = partialcorr(Y);
            R_ij(isnan(R_ij))=0;
            
            R(i,j,1) = R(i,j,1) + R_ij(1,2);
            R(i,j,2) = R(i,j,2) + R_ij(2,3);
            R(i,j,3) = R(i,j,3) + R_ij(1,3);
            R(i,j,4) = R(i,j,4) + R_ij(1,4);
            R(i,j,5) = R(i,j,5) + R_ij(2,4);
            R(i,j,6) = R(i,j,6) + R_ij(3,4);
        end
    end
end

R = R(2:end,2:end,:);
R = R / N_itr;

%% Plots

conns = {'\rho(AB)','\rho(BC)','\rho(CD)'};
names = {'ab','bc','cd'};
idx = [1 2 6];

for i = 1 : 3
   h=figure;
   set(h,'Color','w');
   hh=imagesc(squeeze(R(:,:,idx(i))));
   n = round(length(relC) / 10);
   set(gca,'XTickLabels',relC(n:n:end));
   n = round(length(Ns) / 10);
   set(gca,'YTickLabels',Ns(n+1:n:end));
   set(gca,'FontSize',12);
   hh=title(conns{i});
   set(hh,'FontSize',18);
   hh=xlabel('Relative Noise (C,D)');
   set(hh,'FontSize',16);
   hh=ylabel('Relative Noise (Baseline)');
   set(hh,'FontSize',16);
   caxis([0 1]);
   colorbar;
   
   print(sprintf('par_img_varc_pcorr_%s.png',names{i}),'-dpng','-r0');
end

relbs = [0.5 1 1.5];

% Plot lines
for j = 1 : length(relbs)
    i = find(relC==relbs(j));
    h=figure;
    set(h,'Color','w');
    M = squeeze(R(:,i,idx));
    hh = plot(Ns(2:end),M);
    ylim([0 1]);
    
    set(gca,'FontSize',12);
    
    num = sprintf('%1.1f', relC(i));
    
    hh=title(['w_{c,d} = ' num ' \times baseline']);
    set(hh,'FontSize',18);
    hh=xlabel('Relative Noise (Baseline)');
    set(hh,'FontSize',16);
    hh=ylabel('\rho_{ij}');
    set(hh,'FontSize',16);
    
    hh=legend(conns);
    set(hh,'FontSize',14);
    
    print(sprintf('par_lines_varc_pcorr_wb_%02d.png',relC(i)*10),'-dpng','-r0');
end