%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%           M. Oh, S. Zhao, V. Matveev and F. Nadim (2012)
%    Neuromodulatory changes in short-term synaptic dynamics may be 
%   mediated by two distinct mechanisms of presynaptic calcium entry
%            Journal of Computational Neuroscience. 
%
%         Produces Fig. 2 or Fig. 5 (set Model=1 or Model=2)
%                     Victor Matveev, 2012
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
plfs = 9; tfs = 13; lfs = 12;% FONTSIZES
M = 100; L = 2.8;

filter = 0.5;
Model = 2;
clrs = [ [1 0 0]; 0.9*[1 1 0]; [0 0 1]; [0 0.5 0]; [1 0 1]; [0.5 0.5 0.5]; [1 1 0]; [0 0 0]; [0 1 1]; [0.8 0.8 0.8] ];

minParams = [];
minDev = [];
minCost = [];
normParams = [];
filterParams = [];

if Model == 1
    % tauH tauCa gsyn minfVproc minfSproc tauMmax tauMV tauMS gm hinfV hinfS 
    %   1    2     3     4         5         6      7     8    9  10    11   
    % K4 tauMcntrl minfVcnt minfScnt  gCa 
    % 12    13       14        15     16  
    
    Filter = 8.5;   % Model 1 was run with new cost function, which had larger trace dev contribution
    DataRandGsynModel1;
    minParams = abs(minParams);
    %minParams(:,3) = minParams(:,3) ./ minParams(:,12);
    minParams(:,12) = minParams(:,12) .^ 0.25;
    perm = [12 3 16 9 1 2 13 10 11 14 15 4 5 6 7 8];
    setNum = [1 2];
    myXticks = {'K_{Ca}', 'g_{syn}', 'g_{Ca}', 'g_m', '\tau_h', '\tau_{Ca}', '\tau_m^c', 'V_h', 'S_h', 'V_m', 'S_m',...
                'V_m^p', 'S_m^p', 'A_{\tau}^p', 'V_{\tau}^p','S_{\tau}^p'};
else
    % gMI tauCa gamma kpV kpS gm tauH hinfV hinfS
    %  1   2     3     4   5   6   7   8    9   
    % K4  tauM  minfV minfS gCa  km 
    % 10   11     12    13  14   15
    
    Filter = 6.5;  % Model 2 was run with old cost function, which had smaller trace dev contribution
    DataRandGsynModel2;
    minParams = abs(minParams);
    %minParams(:,3) = minParams(:,3) ./ minParams(:,10);
    minParams(:,10) = minParams(:,10) .^ 0.25;
    minParams = abs(minParams);
    perm = [10 3 14 6 1 7 2 11 8 9 12 13 4 5];
    setNum = [1 8] ;%[18 12];
    myXticks = {'K_{Ca}', 'g_{syn}', 'g_{Ca}', 'g_m', 'g_{MI}', '\tau_h', '\tau_{Ca}', '\tau_{m}', 'V_{h}', 'S_{h}', ...
                'V_{m}', 'S_{m}', 'V^{+}_{}', 'S^{+}_{}'};
end;


SZ = size(minParams);
meanParams = mean(minParams);
N = SZ(1);
nPars = length(perm);

for k = 1 : SZ(1)
    if minCost(k) < Filter && ~sum(minDev(k,2:5) > 0.2)
        filterParams = [filterParams; minParams(k, perm) ];
    end
end
for k = 1 : size(filterParams,1)
        normParams = [normParams; filterParams(k,:) ./ mean(filterParams) ];
end

CVarray = sqrt( var(filterParams) ) ./ mean(filterParams);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(Model);
subplot(121);
%set(gca,'fontsize',12);
hold on;

width=6;
chart = zeros( nPars+1, M );

for ii = 1 : size(normParams,1)
    for jj = 1 : nPars 
        yy = normParams(ii, jj) / L * M;
        for kk = 1 : M 
            chart(jj, kk) = chart(jj, kk) + exp(-(kk-yy)^2/width);;
       end
    end
end

maxChart = max(max(chart));
chart =  maxChart - chart;
aveChart = mean( mean(chart) );

pcolor( chart' );
shading flat;

colormap gray;
caxis([0.6*aveChart maxChart]);

hold on;

for ii = 1 : nPars
    plot([ii ii], [0 M], '-', 'color', [1 1 1], 'linewidth', 2);
    text(ii+0.15, -2, myXticks{ii}, 'fontsize', plfs);
end

set(gca,'Xtick', []);
set(gca,'Ytick', []);
%set(gca,'box','off');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

plot([1.05 1.05], [0 M], 'k-', 'linewidth', 2);
plot([0 nPars+1], 0.5+[1 1]*M/L, 'b--', 'linewidth', 2);

[Y, I] = sort(minCost);
ii = 1;
for s = setNum  %length(Sets)
    plot( [1 : nPars] + 0.5, minParams( I(s), perm) ./ mean(filterParams) * M / L , ...
          'o-', 'linewidth', 2, 'color', clrs(ii,:), 'markerfacecolor', clrs(ii,:), 'markeredgecolor', clrs(ii,:) );
    ii = ii + 1;
end

ylabel('Normalized value', 'fontsize', lfs);
title('{\bf (A)}   Normalized parameter values at minimum', 'fontsize', tfs);
axis([1, nPars+1, 1, M]);

for ii = 1 : nPars
    if CVarray(ii) < 0.4
        text(ii+0.5, 0.84*M, sprintf('CV = %.2g', CVarray(ii) ), 'color', [0 0 1], 'rotation', 90);
    else
        text(ii+0.5, 0.84*M, sprintf('CV = %.2g', CVarray(ii) ), 'rotation', 90);
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

CORR = zeros(nPars, nPars);

for i1 = 1 : nPars
    for i2 = i1 : nPars
        CORR(i1,i2) = (normParams(:,i1) - 1)' * (normParams(:,i2) - 1) / ...
                  sqrt(sum((normParams(:,i1) - 1).^2) * sum((normParams(:,i2) - 1).^2) );
    end
end

CORR = CORRCOEF(filterParams);
CORR = CORR .* (abs(CORR) > filter); 
%for ii = 1 : nPars - 1
%    for jj = ii + 1 : nPars
%        CORR(jj,ii) = 1;
%    end
%end;

chain = cell(nPars,1);
chainNum = 0;
fprintf('\n\n Correlation structure: \n\n');

for i1 = 1 : (nPars - 1)
    for i2 = i1+1 : nPars
        if CORR(i1,i2) > filter
            fprintf('Positive correlation pair: %d %d\n', i1, i2);
            match = 0;
            for id = 1 : chainNum
                if sum( chain{id} == i1 ) && sum(chain{id} == i2 )
                    match = 1;
                    break;
                elseif sum( chain{id} == i1 ) 
                    chain{id} = [ chain{id}, i2 ];
                    match = 1;
                    break;
                elseif sum( chain{id} == i2 ) 
                    chain{id} = [ chain{id}, i1 ];
                    match = 1;
                    break;
                end
            end
            if ~match
                chainNum = chainNum + 1;
                chain{chainNum} = [i1 i2];
            end
        end
    end
end

for i = 1 : chainNum
    fprintf('Correlation chain #%d', i);
    disp( sort(chain{i}) );
end;

for i1 = 1 : (nPars - 1)
    for i2 = i1 + 1 : nPars
        if CORR(i1,i2) < -filter
            fprintf('Negative correlation pair: %d %d\n', i1, i2);
        end
    end
end

CV = zeros(nPars,1);
for ii = 1 : nPars
    CV(ii) = sqrt(sum((normParams(:, ii) - 1).^2));
end

subplot(122);
%set(gca,'fontsize',12);
PCORR = [CORR, zeros(nPars,1) ];
PCORR = [PCORR; zeros(1,nPars+1) ];
pcolor(1:(nPars+1), 1:(nPars+1), PCORR);
colormap gray;
colorbar('NorthOutside');
hold on;

shading flat;
set(gca, 'XTick', []);
set(gca, 'YTick', []);

for ii = 1 : nPars
    text(ii+0.22, 0.2, myXticks{ii}, 'fontsize', plfs);
    text(nPars+1.2, ii+0.4, myXticks{ii}, 'fontsize', plfs);
    plot([ii ii], [0 nPars+1], '-', 'color', [1 1 1], 'linewidth', 1);
    plot([0 nPars+1], [ii ii],  '-', 'color', [1 1 1], 'linewidth', 1);
end

title('{\bf (B)}   Parameter Correlations', 'fontsize', tfs);
fprintf('%.3g ',minParams(I(setNum(1)),perm) );

