% Written by
%   Alessandro Motta <alessandro.motta@brain.mpg.de>
clear;

%% Configuration
rnnConfigs = struct;

rnnConfigs(1).tag = 'DISC';
rnnConfigs(1).rootDir = fullfile( ...
    '/tmpscratch', 'amotta', 'klinger', ...
    '2019-07-23-inh-exc-relu-rnn-with-pruning');
rnnConfigs(1).measureCsvDir = fullfile( ...
    rnnConfigs(1).rootDir, 'measures', '20200415T170632');
rnnConfigs(1).perfCsvDir = fullfile( ...
    rnnConfigs(1).rootDir, 'exported-csvs', '20190926T160505');
rnnConfigs(1).connectomeDir = fullfile( ...
    rnnConfigs(1).rootDir, 'connectomes', '20200422T114417');

rnnConfigs(2).tag = 'MEM';
rnnConfigs(2).rootDir = fullfile( ...
    '/tmpscratch', 'amotta', 'klinger', ...
    '2020-03-01-inh-exc-relu-memory-rnn-with-pruning');
rnnConfigs(2).measureCsvDir = fullfile( ...
    rnnConfigs(2).rootDir, 'measures', '20200415T142543');
rnnConfigs(2).perfCsvDir = fullfile( ...
    rnnConfigs(2).rootDir, 'tfsummaries', '20200318T175925');
rnnConfigs(2).connectomeDir = fullfile( ...
    rnnConfigs(2).rootDir, 'connectomes', '20200422T101227');

curRoot  = fullfile( ...
    '/mnt', 'mpibr', 'data', 'Personal', 'mottaa', 'Klinger');
modelSummaryXlsFile = fullfile(curRoot, ...
    '2019-10-28-Measure-Distributions', 'description_noisefree_v1.xlsx');
modelRawXlsDir = fullfile(curRoot, ...
    '2020-01-13-Raw-Measure-Data-Points', 'samples_noisefree');

vars = struct;
vars(1).rnnName = 'relative_reciprocity_ee';
vars(1).modelName = 'eeRelativeReciprocity';
vars(1).yLabel = 'r_{ee}';
vars(1).yTicks = 10 .^ (0:2);
vars(1).yLim = 10 .^ [-0.5, 2.5];
vars(1).yScale = 'log';

vars(2).rnnName = 'relative_reciprocity_ei';
vars(2).modelName = 'eiRelativeReciprocity';
vars(2).yLabel = 'rr_{ei}';
vars(2).yTicks = 10 .^ (0:2);
vars(2).yLim = 10 .^ [-0.5, 2.5];
vars(2).yScale = 'log';

vars(3).rnnName = 'relative_reciprocity_ie';
vars(3).modelName = 'ieRelativeReciprocity';
vars(3).yLabel = 'rr_{ie}';
vars(3).yTicks = 10 .^ (0:2);
vars(3).yLim = 10 .^ [-0.5, 2.5];
vars(3).yScale = 'log';

vars(4).rnnName = 'relative_reciprocity_ii';
vars(4).modelName = 'iiRelativeReciprocity';
vars(4).yLabel = 'rr_{ii}';
vars(4).yTicks = 10 .^ (0:2);
vars(4).yLim = 10 .^ [-0.5, 2.5];
vars(4).yScale = 'log';

vars(5).rnnName = 'relative_recurrency_ee';
vars(5).modelName = 'relativeCycles5';
vars(5).yLabel = 'r^{(5)}';
vars(5).yTicks = 10 .^ (0:5);
vars(5).yLim = 10 .^ [-0.5, 5];
vars(5).yScale = 'log';

vars(6).rnnName = 'in_out_degree_correlation_ee';
vars(6).modelName = 'excInOutDegreeCorrelation';
vars(6).yLabel = 'r_{i/o}';
vars(6).yTicks = (-1):1;
vars(6).yLim = [-1, +1];
vars(6).yScale = 'linear';

% See "circuit constraints" section of methods
numExc = 1800; connExc = 0.2;
numInh = 200; connInh = 0.6;
numAll = numExc + numInh;

numConns = ...
    numExc * (connExc * (numAll - 1)) ...
  + numInh * (connInh * (numAll - 1));
fracConn = ...
    numConns / (numAll * (numAll - 1));

clear cur*;
info = Util.runInfo();
Util.showRunInfo(info);

%% Loading network measure data
clear cur*;

rnnTypes = categorical(cell(0, 1));
rnnPerfs = cell(0, 1);
rnnConnectomes = cell(0, 2);
rnnMeasuresThresh = cell(0, 1);
rnnMeasures = cell(0, 1);

for curIdx = 1:numel(rnnConfigs)
    curRnnConfig = rnnConfigs(curIdx);
    
    curNames = dir(fullfile(curRnnConfig.perfCsvDir, '*.csv'));
    curNames = reshape({curNames.name}, [], 1);
    
    curMeasureCsvFiles = fullfile(curRnnConfig.measureCsvDir, curNames);
    curPerfCsvFiles = fullfile(curRnnConfig.perfCsvDir, curNames);
    
    switch curRnnConfig.tag
        case 'DISC'
            curPerfMetrics = {'loss', 'accuracy'};
            curColors = [0, 0, 0; 1, 1, 1];
        case 'MEM'
            curPerfMetrics = {'loss'};
            curColors = [0, 0, 0; 1, 0, 1];
        otherwise
            error('Invalid tag "%s"', curRnnConfig.tag);
    end
    
    curInterp = linspace(0, 1, numel(curPerfCsvFiles) + 2);
    curInterp = reshape(curInterp(2:(end - 1)), [], 1);
    
    curColors = ( ...
        (curColors(1, :) .^ 2.2) .* (0 + curInterp) ...
      + (curColors(2, :) .^ 2.2) .* (1 - curInterp)) .^ (1 / 2.2);
    rnnConfigs(curIdx).colors = curColors;
    
    curTypes = categorical(repelem( ...
        {curRnnConfig.tag}, numel(curPerfCsvFiles)));
   [~, curPerfs] = cellfun( ...
        @(p) RNN.loadPerformance(p, curPerfMetrics), ...
        curPerfCsvFiles, 'UniformOutput', false);
    
    curConnectomes = cell(numel(curMeasureCsvFiles), 2);
    curMeasuresThresh = cell(size(curMeasureCsvFiles));
    curMeasures = cell(size(curMeasureCsvFiles));
    for curRnnIdx = 1:numel(curMeasureCsvFiles)
        curRnnName = curNames{curRnnIdx};
       [~, curRnnName] = fileparts(curRnnName);
       
        curRnnMeasures = readtable( ...
            curMeasureCsvFiles{curRnnIdx}, ...
            'Delimiter', ',', 'ReadVariableNames', true);
        curMeasuresThresh{curRnnIdx} = curRnnMeasures;
        
        % NOTE(amotta): Let's drop the measures obtained after thresholding
        % the weighted connectome at a particular percentile.
        curMask = curRnnMeasures.threshold_percentile == 0;
        curRnnMeasures = curRnnMeasures(curMask, :);
        curMeasures{curRnnIdx} = curRnnMeasures;
        
        % NOTE(amotta): Find the checkpoint with the overall pairwise
        % connectivity that is closest to the experimental data and then
        % load its binary binary and weighted connectomes.
        curCkptIdx = curRnnMeasures.connectivity_aa;
       [~, curCkptIdx] = min(abs(curCkptIdx - fracConn));
       
        curCkptName = curRnnMeasures.ckpt_file{curCkptIdx};
        curCkptName = strsplit(curCkptName, filesep);
        curCkptName = curCkptName{end};
        
        curConnFile = fullfile( ...
            curRnnConfig.connectomeDir, ...
            curRnnName, sprintf('%s.hdf5', curCkptName));
        curConnMask = logical(h5read(curConnFile, '/mask'));
        curConnWeights = h5read(curConnFile, '/weights');
        
        curConnectomes{curRnnIdx, 1} = curConnMask;
        curConnectomes{curRnnIdx, 2} = curConnWeights;
    end
    
    rnnTypes = [rnnTypes; curTypes(:)]; %#ok
    rnnPerfs = [rnnPerfs; curPerfs(:)]; %#ok
    rnnConnectomes = [rnnConnectomes; curConnectomes]; %#ok
    rnnMeasuresThresh = [rnnMeasuresThresh; curMeasuresThresh(:)]; %#ok
    rnnMeasures = [rnnMeasures; curMeasures(:)]; %#ok
end

[modelSummaryT, modelRawData] = ...
    RNN.loadModelMeasures(modelSummaryXlsFile, modelRawXlsDir);
modelPlotConfigT = RNN.buildModelPlotConfig();

%% Plot performance versus connectivity
clear cur*;
curAxesY = {'left', 'right'};

% NOTE(amotta): This is the connectivity range in which we consider all
% RNNs to have learned (and not yet "forgotten") the task. It was
% determined manually by inspection of the figure.
curEvalConnRange = flip([0.8, 0.00375]);

curFig = figure();
curAx = axes(curFig);

for curConfigIdx = flip(1:numel(rnnConfigs))
    curConfig = rnnConfigs(curConfigIdx);
    
    switch curConfig.tag
        case 'DISC'
            curMetric = 'accuracy';
            curSummaryFun = @max;
            curDirY = 'normal';
        case 'MEM'
            curMetric = 'loss';
            curSummaryFun = @min;
            curDirY = 'reverse';
    end
    
    curPerfs = rnnTypes == curConfig.tag;
    curPerfs = rnnPerfs(curPerfs);
    curColors = curConfig.colors;
    
    curAxisY = curAxesY{curConfigIdx};
    yyaxis(curAx, curAxisY);
    hold(curAx, 'on');
    
    for curRnnIdx = 1:numel(curPerfs)
        curPerfT = curPerfs{curRnnIdx};
        curColor = curColors(curRnnIdx, :);
        
       [curX, ~, curY] = unique(curPerfT.weightCount);
        curY = accumarray(curY, curPerfT.(curMetric), [], curSummaryFun);
        curX = curX / (numAll * (numAll - 1));
        
        curMeanPerf = ...
            (curEvalConnRange(1) < curX) ...
          & (curX < curEvalConnRange(2));
        curMeanPerf = mean(curY(curMeanPerf));
        
        fprintf( ...
            '%s RNN %d. Average performance: %f\n', ...
            curConfig.tag, curRnnIdx, curMeanPerf);

        plot(curAx, curX, curY, '-', 'Color', curColor, 'LineWidth', 2);
    end
    
    ylim(curAx, [0, 1]);
    yticks(curAx, 0:0.2:1);
    curAx.YTickLabel(2:(end - 1)) = {''};
    
    ylabel(curAx, curMetric);
    curAx.YDir = curDirY;
end

axis(curAx, 'square');

curAx.XScale = 'log';
curAx.XDir = 'reverse';

xlabel(curAx, 'Connectivity');
xlim(curAx, flip([1, 1E-5]));
xticks(curAx, 10 .^ (-5:0));
xticklabels(curAx, {'', '', '0.001', '', '', '1'});

Figure.config(curFig, info);
curFig.Position(3:4) = [300, 190];

%% Plot network measures versus connectivity
clear cur*;

curDataLim = flip([1E0, 10^(-3.5)]);
curLimX = curDataLim;

curTicksX = 10 .^ (-3:0);
curTickLabelsX = arrayfun( ...
    @(x) strrep(num2str(x, '%g'), '0.', '.'), ...
    curTicksX, 'UniformOutput', false);

curColors = cat(1, rnnConfigs.colors);
curFig = figure();

for curVarIdx = 1:numel(vars)
    curVar = vars(curVarIdx);
    curModelVarName = curVar.modelName;
    
    curAx = subplot(1, numel(vars), curVarIdx);
    hold(curAx, 'on');
    
    % Indicate connectivity fractions consistent with mouse barrel
    plot(curAx, repelem(fracConn, 2), curVar.yLim, 'k--');
    
    for curRnnIdx = 1:numel(rnnMeasures)
        curX = rnnMeasures{curRnnIdx}.connectivity_aa;
        curY = rnnMeasures{curRnnIdx}.(curVar.rnnName);
        curColor = curColors(curRnnIdx, :);
        plot(curAx, curX, curY, 'Color', curColor, 'LineWidth', 2);
    end
end

curFig.Position(3:4) = [850, 190];
curFig.Color = 'white';

curAxes = findobj(curFig, 'Type', 'Axes');
curAxes = flip(curAxes);

arrayfun(@(a) xticks(a, curTicksX), curAxes);
arrayfun(@(a) xticklabels(a, curTickLabelsX), curAxes);
arrayfun(@(a, c) set(a, 'YScale', c.yScale), curAxes(:), vars(:));
arrayfun(@(a, c) ylim(a, c.yLim), curAxes(:), vars(:));
arrayfun(@(a, c) yticks(a, c.yTicks), curAxes(:), vars(:));
arrayfun(@(a, c) ylabel(a, c.yLabel), curAxes(:), vars(:));

set(curAxes, ...
    'TickDir', 'out', ...
    'XDir', 'reverse', ...
    'XScale', 'log', ...
   {'XLim'}, {curLimX}, ...
   {'TickLength'}, {[0.05, 0.025]}, ...
   {'XMinorTick'}, {'on'}, ...
   {'YMinorTick'}, {'on'});

arrayfun(@(a) axis(a, 'square'), curAxes);
xlabel(curAxes(1), 'Connectivity');

annotation( ...
    curFig, ...
    'textbox', [0, 0.9, 1, 0.1], ...
	'String', {info.filename; info.git_repos{1}.hash}, ...
    'EdgeColor', 'none', 'HorizontalAlignment', 'center');

%% Plot raw data points
clear cur*;

curVars = vars;
[curVars.yScale] = deal('linear');

% NOTE(amotta): From panel 3c
curVars(1).yLim = [0, 4]; curVars(1).yTicks = 0:0.1:4;
curVars(2).yLim = [0, 2]; curVars(2).yTicks = 0:0.1:2;
curVars(3).yLim = [0, 2]; curVars(3).yTicks = 0:0.1:2;
curVars(4).yLim = [0, 2]; curVars(4).yTicks = 0:0.1:2;
curVars(5).yLim = [0, 10]; curVars(5).yTicks = 0:1:10;
curVars(6).yLim = [-1, +1]; curVars(6).yTicks = (-1):0.1:1;

curRnnDataT = nan(numel(rnnMeasures), numel(curVars));

for curRnnIdx = 1:numel(rnnMeasures)
    curX = rnnMeasures{curRnnIdx}.connectivity_aa;
   [~, curIdx] = min(abs(curX - fracConn));
    
    for curVarIdx = 1:numel(curVars)
        curVarName = curVars(curVarIdx).rnnName;
        curDatum = rnnMeasures{curRnnIdx}.(curVarName);
        curDatum = curDatum(curIdx);
        
        curRnnDataT(curRnnIdx, curVarIdx) = curDatum;
    end
end

curRnnDataT = array2table(curRnnDataT, ...
    'VariableNames', {curVars.modelName});

curRawDataT = modelRawData;
curRawDataT.discrnn = curRnnDataT(rnnTypes == 'DISC', :);
curRawDataT.memrnn = curRnnDataT(rnnTypes == 'MEM', :);

curModelT = [modelPlotConfigT; { ...
    'DISCRNN', 'DISC-RNN', [0, 0, 0]; ...
    'MEMRNN',  'MEM-RNN',  [1, 0, 1]}];

dup = @(x) [x, x];
curFig = figure();
for curVarIdx = 1:numel(curVars)
    curVar = curVars(curVarIdx);
    curModelVarName = curVar.modelName;
    
    curAx = subplot(1, numel(curVars), curVarIdx);
    hold(curAx, 'on');
    
    for curModelIdx = 1:height(curModelT)
        curModelName = lower(curModelT.tag{curModelIdx});
        curModelSummaryT = curRawDataT.(curModelName);
        curColor = curModelT.color(curModelIdx, :);
        
        curY = curModelSummaryT.(curModelVarName);
        curX = curModelIdx + 0.5 * (rand(numel(curY), 1) - 0.5);
        scatter(curAx, curX, curY, 24, curColor, '.');
        
        curX = curModelIdx + 0.75 * [-0.5, +0.5];
        curBox = prctile(curY, [0, 25, 50, 75, 100]);
        curWhiskerX = curModelIdx + 0.75 * [-0.5, +0.5];
        
        rectangle(curAx, ...
            'Position', [ ...
                curX(1), curBox(2), ...
                curX(end) - curX(1), ...
                curBox(4) - curBox(2)]);
        plot(curAx, ...
            curX, dup(curBox(3)), ...
            'Color', 'black');
        
        plot(curAx, curWhiskerX, dup(curBox(1)), 'Color', 'black');
        plot(curAx, curWhiskerX, dup(curBox(end)), 'Color', 'black');
        
        plot( ...
            curAx, dup(curModelIdx), [curBox(4), curBox(end)], ...
            'Color', 'black', 'LineStyle', '--');
        plot( ...
            curAx, dup(curModelIdx), [curBox(1), curBox(2)], ...
            'Color', 'black', 'LineStyle', '--');
    end
end

curFig.Position(3:4) = [850, 190];
curFig.Color = 'white';

curAxes = findobj(curFig, 'Type', 'Axes');
curAxes = flip(curAxes);

curLimX = [0.5, height(curModelT) + 0.5];
curTicksX = 1:height(curModelT);

arrayfun(@(a, c) set(a, 'YScale', c.yScale), curAxes(:), curVars(:));
arrayfun(@(a, c) ylim(a, c.yLim), curAxes(:), curVars(:));
% arrayfun(@(a, c) yticks(a, c.yTicks), curAxes(:), curConfigs(:));
arrayfun(@(a, c) ylabel(a, c.yLabel), curAxes(:), curVars(:));

set(curAxes, ...
    'TickDir', 'out', ...
   {'XLim'}, {curLimX}, ...
   {'XTick'}, {curTicksX}, ...
   {'TickLength'}, {[0.05, 0.025]}, ...
   {'YMinorTick'}, {'on'});

arrayfun(@(a) axis(a, 'square'), curAxes);
xlabel(curAxes(1), 'Model');

annotation( ...
    curFig, ...
    'textbox', [0, 0.9, 1, 0.1], ...
	'String', {info.filename; info.git_repos{1}.hash}, ...
    'EdgeColor', 'none', 'HorizontalAlignment', 'center');

%% Quantitative evaluation
curVarNames = {curVars.modelName};

curEvalT = curModelT(:, {'tag'});
curEvalT.Properties.RowNames = curModelT.title;

curEvalT.data = cellfun( ...
    @(t) curRawDataT.(lower(t)), ...
    curEvalT.tag, 'UniformOutput', false);

[curMeanT, curVarT] = cellfun( ...
    @(t) deal( ...
        mean(t{:, curVarNames}, 1), ...
        var(t{:, curVarNames}, 1)), ...
    curEvalT.data, 'UniformOutput', false);

curMeanT = array2table(cell2mat(curMeanT), ...
    'VariableNames', curVarNames, 'RowNames', curModelT.title);
curVarT = array2table(cell2mat(curVarT), ...
    'VariableNames', curVarNames, 'RowNames', curModelT.title);
curStdT = curVarT; curStdT{:, :} = sqrt(curStdT{:, :});

% Statistical evaluation
curModelNames = {'DISC-RNN', 'MEM-RNN'};
curPvals = nan(1, numel(curVarNames));
for curVarIdx = 1:numel(curVarNames)
    curVarName = curVarNames{curVarIdx};
    
    curMeans = curMeanT{curModelNames, curVarName};
    curStds = curStdT{curModelNames, curVarName};
    
    curData = cellfun( ...
        @(t) t{:, curVarName}, ...
        curEvalT.data(curModelNames), ...
        'UniformOutput', false);
   [~, curPval] = kstest2(curData{:});
    
    fprintf('%s: %.2f+-%.2f vs. %.2f+-%.2f, p=%.3f (%s vs. %s)\n', ...
        curVarName, curMeans(1), curStds(1), ...
        curMeans(2), curStds(2), curPval, ...
        curModelNames{:});
end

%% Reproducing figure 4.4 of Emmanuel's PhD thesis
clear cur*;
rng(2);

curModelNames = modelSummaryT.Properties.RowNames;

% Collecting data
curRnnVarNames = {vars.modelName};
modelTsneT = cell(height(modelSummaryT), 1);
for curModelIdx = 1:height(modelSummaryT)
    curModelName = curModelNames{curModelIdx};
    
    curModelSummaryT = lower(curModelName);
    curModelSummaryT = modelRawData.(curModelSummaryT);
    curModelSummaryT = curModelSummaryT(:, curRnnVarNames);
    
    curModelSummaryT.model(:) = categorical({curModelName});
    curModelSummaryT.connectivity(:) = fracConn;
    curModelSummaryT = circshift(curModelSummaryT, 2, 2);
    
    modelTsneT{curModelIdx} = curModelSummaryT;
end

% tSNE
curTsnes = cat(1, modelTsneT{:});
curTsnes = tsne(curTsnes{:, 3:end});
curTsnes = mat2cell(curTsnes, cellfun(@height, modelTsneT), 2);

% Visualization using tSNE
curFig = figure();
curAx = axes(curFig);
hold(curAx, 'on');

for curModelIdx = 1:height(modelSummaryT)
    curColor = curModelNames{curModelIdx};
    curColor = strcmpi(modelPlotConfigT.tag, curColor);
    curColor = modelPlotConfigT.color(curColor, :);
    
    curTsne = curTsnes{curModelIdx};
    
    scatter(curAx, ...
        curTsne(:, 1), curTsne(:, 2), ...
        [], curColor, 'filled');
end

xlabel(curAx, 'First t-SNE dimension');
ylabel(curAx, 'Second t-SNE dimension');

axis(curAx, 'square');
curLims = [xlim(curAx), ylim(curAx)];
curLims = [-1, +1] * max(abs(curLims));
xlim(curAx, curLims);
ylim(curAx, curLims);

[~, curLeg] = ismember(curModelNames, modelPlotConfigT.tag);
curLeg = modelPlotConfigT.title(curLeg);

legend(curAx, curLeg, 'Location', 'EastOutside');

Figure.config(curFig, info);
curFig.Position(3:4) = [420, 260];

%% Visualize training of RNNs
clear cur*;
rng(8);

curRnnVarNames = [{'connectivity_aa'}, vars.rnnName];
curModelVarNames = [{'connectivity'}, vars.modelName];

rnnTsneT = cell(numel(rnnMeasures), 1);
for curRnnIdx = 1:numel(rnnMeasures)
    curRnnType = rnnTypes(curRnnIdx);
    
    curRnnName = rnnTypes(1:curRnnIdx);
    curRnnName = sum(curRnnName == curRnnType);
    curRnnName = sprintf('RNN-%s-%d', curRnnType, curRnnName);
    
    curRnnDataT = rnnMeasures{curRnnIdx};
    curRnnDataT = sortrows(curRnnDataT, 'connectivity_aa', 'descend');
    
    % NOTE(amotta); Only show up to biologicall plausible connectivity.
    curRnnDataT = curRnnDataT(curRnnDataT.connectivity_aa >= fracConn, :);
    
    % NOTE(amotta): In-ou-degree correlation is NaN if fully connected
    curRnnDataT(1, :) = [];
    
    curRnnDataT = curRnnDataT(:, curRnnVarNames);
    curRnnDataT.Properties.VariableNames = curModelVarNames;
    
    curRnnDataT.model(:) = categorical({curRnnName});
    curRnnDataT = circshift(curRnnDataT, 1, 2);
    rnnTsneT{curRnnIdx} = curRnnDataT;
end

curRnnNames = cellfun( ...
    @(t) char(t.model(1)), rnnTsneT, 'UniformOutput', false);
curRnnColors = reshape({rnnConfigs.colors}, [], 1);
curRnnColors = num2cell(cell2mat(curRnnColors), 2);

curTsneT = [modelTsneT; rnnTsneT];
curPlotConfigT = [modelPlotConfigT; ...
    curRnnNames, curRnnNames, curRnnColors];

curNames = cellfun( ...
    @(data) char(data.model(1)), ...
    curTsneT, 'UniformOutput', false);
curNames = reshape(curNames, 1, []);

% tSNE
curTsnes = cat(1, curTsneT{:});
curTsnes = tsne(curTsnes{:, 3:end});
curTsnes = mat2cell(curTsnes, cellfun(@height, curTsneT), 2);

% Visualization using tSNE
curFig = figure();
curAx = axes(curFig);
hold(curAx, 'on');

for curDataIdx = 1:numel(curTsneT)
    curName = char(curTsneT{curDataIdx}.model(1));
    
    curColor = strcmpi(curPlotConfigT.tag, curName);
    curColor = curPlotConfigT.color(curColor, :);
    
    curTsne = curTsnes{curDataIdx};
    
    if startsWith(curName, 'RNN')
        curSizes = repelem(12, size(curTsne, 1));
        curSizes(end) = 48;
        
        plot(curAx, ...
            curTsne(:, 1), curTsne(:, 2), ...
            'Color', curColor, 'LineWidth', 1.5);
        scatter(curAx, ...
            curTsne(:, 1), curTsne(:, 2), ...
            curSizes, curColor, 'filled');
    else
        scatter(curAx, ...
            curTsne(:, 1), curTsne(:, 2), ...
            12, curColor, 'filled');
    end
end

curIsLine = arrayfun(@(p) isgraphics(p, 'Line'), curAx.Children);
[curIsLine, curSortIds] = sort(curIsLine, 'ascend');
curAx.Children = curAx.Children(curSortIds);

xlabel(curAx, 'First t-SNE dimension');
ylabel(curAx, 'Second t-SNE dimension');
axis(curAx, 'equal');

[~, curLeg] = ismember(curNames, curPlotConfigT.tag);
curLeg = curPlotConfigT.title(curLeg);

legend( ...
    flip(curAx.Children(not(curIsLine))), ...
    curLeg, 'Location', 'EastOutside');

Figure.config(curFig, info);
curFig.Position(3:4) = [950, 620];

%% Compare RNNs in t-SNE space over large connectivity range
clear cur*;
rng(9);

curRnnVarNames = [{'connectivity_aa'}, vars.rnnName];
curModelVarNames = [{'connectivity'}, vars.modelName];

% NOTE(amotta): Restrict to range covered by all RNNs
curMinConn = max(cellfun(@(m) min(m.connectivity_aa), rnnMeasures));

curTsneT = cell(numel(rnnMeasures), 1);
for curRnnIdx = 1:numel(rnnMeasures)
    curRnnType = rnnTypes(curRnnIdx);
    
    curRnnName = rnnTypes(1:curRnnIdx);
    curRnnName = sum(curRnnName == curRnnType);
    curRnnName = sprintf('RNN-%s-%d', curRnnType, curRnnName);
    
    curRnnDataT = rnnMeasures{curRnnIdx};
    curRnnDataT = sortrows(curRnnDataT, 'connectivity_aa', 'descend');
    
    % NOTE(amotta); Only show up to biologicall plausible connectivity.
    curRnnDataT = curRnnDataT(curRnnDataT.connectivity_aa >= curMinConn, :);
    
    % NOTE(amotta): In-ou-degree correlation is NaN if fully connected
    curRnnDataT(1, :) = [];
    
    curRnnDataT = curRnnDataT(:, curRnnVarNames);
    curRnnDataT.Properties.VariableNames = curModelVarNames;
    
    curRnnDataT.model(:) = categorical({curRnnName});
    curRnnDataT = circshift(curRnnDataT, 1, 2);
    curTsneT{curRnnIdx} = curRnnDataT;
end

curRnnNames = cellfun( ...
    @(t) char(t.model(1)), rnnTsneT, 'UniformOutput', false);
curRnnColors = reshape({rnnConfigs.colors}, [], 1);
curRnnColors = num2cell(cell2mat(curRnnColors), 2);

curPlotConfigT = [modelPlotConfigT([], :); ...
    curRnnNames, curRnnNames, curRnnColors];

curNames = cellfun( ...
    @(data) char(data.model(1)), ...
    curTsneT, 'UniformOutput', false);
curNames = reshape(curNames, 1, []);

% tSNE
curTsnes = cat(1, curTsneT{:});
curTsnes = tsne(curTsnes{:, 3:end});
curTsnes = mat2cell(curTsnes, cellfun(@height, curTsneT), 2);

% Visualization using tSNE
curFig = figure();
curAx = axes(curFig);
hold(curAx, 'on');

for curDataIdx = 1:numel(curTsneT)
    curTsne = curTsnes{curDataIdx};
    curName = char(curTsneT{curDataIdx}.model(1));
    
    curColor = strcmpi(curPlotConfigT.tag, curName);
    curColor = curPlotConfigT.color(curColor, :);
    
    curSizes = repelem(12, size(curTsne, 1));
    curSizes(end) = 48;

    plot(curAx, ...
        curTsne(:, 1), curTsne(:, 2), ...
        'Color', curColor, 'LineWidth', 1.5);
    scatter(curAx, ...
        curTsne(:, 1), curTsne(:, 2), ...
        curSizes, curColor, 'filled');
end

curIsLine = arrayfun(@(p) isgraphics(p, 'Line'), curAx.Children);
[curIsLine, curSortIds] = sort(curIsLine, 'ascend');
curAx.Children = curAx.Children(curSortIds);

xlabel(curAx, 'First t-SNE dimension');
ylabel(curAx, 'Second t-SNE dimension');
axis(curAx, 'equal');

[~, curLeg] = ismember(curNames, curPlotConfigT.tag);
curLeg = curPlotConfigT.title(curLeg);

legend( ...
    flip(curAx.Children(not(curIsLine))), ...
    curLeg, 'Location', 'EastOutside');

Figure.config(curFig, info);
curFig.Position(3:4) = [950, 620];

%% Quantitative evaluation of separability
curMaxConn = min(cellfun(@(t) t.connectivity(1), curTsneT));
curMinConn = max(cellfun(@(t) t.connectivity(end), curTsneT));

curConns = cat(1, curTsneT{:});
curConns = sort(curConns.connectivity, 'descend');
curConns = curConns(curMinConn <= curConns & curConns <= curMaxConn);

curConnCount = numel(curConns);
curDPrimes = nan(curConnCount, 1);

[~, ~, curTypeIds] = unique(rnnTypes);
curIntraMask = reshape(curTypeIds, [], 1) == reshape(curTypeIds, 1, []);

for curConnIdx = 1:curConnCount
    curConn = curConns(curConnIdx);
    
    curConnData = cat(1, curTsneT{:});
    curConnData.absDiff = abs(curConnData.connectivity - curConn);
    curConnData = sortrows(curConnData, 'absDiff');
    
   [curModels, curRows] = unique(curConnData.model);
    curConnData = curConnData(curRows, :);
    curConnData = curConnData{:, 3:end};
    
    curNormConsts = prctile(curConnData, [20, 80], 1);
    curNormConsts = diff(curNormConsts, 1, 1);
    
    curDists = @(a, b) sum(abs(a - b) ./ curNormConsts, 2);
    curDists = squareform(pdist(curConnData, curDists));
    
    curUniDists = unique(curDists);
    curAccuracies = nan(size(curUniDists));
    
    % NOTE(amotta): The d' ("dee prime") sensitivity index can be
    % approximated by sqrt(2) * Z-transform(area under receiver curve).
    % See https://en.wikipedia.org/w/index.php?title=Sensitivity_index&oldid=963482225
   [~, ~, ~, curAuc] = perfcurve(curIntraMask(:), -curDists(:), true);
    curAucDPrime = sqrt(2) * norminv(curAuc);
    curDPrimes(curConnIdx) = curAucDPrime;
end

curFig = figure();
curAx = axes(curFig);
axis(curAx, 'square');
hold(curAx, 'on');

plot(curAx, ...
    curConns, curDPrimes, ...
    'Color', 'black', 'LineWidth', 2);

xlabel(curAx, 'Connectivity');
curAx.XDir = 'reverse';
curAx.XScale = 'log';

ylim(curAx, [0, 2]);
ylabel(curAx, 'Separation of RNNs');

curFig.Position(3:4) = 330;
Figure.config(curFig, info);

curMaxConn = 0.11;
curDPrimeMean = mean(curDPrimes(curConns <= curMaxConn));
curDPrimeStd = std(curDPrimes(curConns <= curMaxConn));

fprintf( ...
    'd'' for p below %.2f is %.2f+-%.2f\n', ...
    curMaxConn, curDPrimeMean, curDPrimeStd);

%% Plot histogram of connection strengths
clear cur*;

curRnnIds = [2, 5];
curBinEdges = linspace(-2.5, -0, 51);

curFig = figure();
curAx = axes(curFig);
axis(curAx, 'square');
hold(curAx, 'on');

for curRnnIdx = 1:numel(curRnnIds)
    curRnnId = curRnnIds(curRnnIdx);
   [curMask, curWeights] = deal(rnnConnectomes{curRnnId, :});
    assert(not(any(curWeights(not(curMask)))));
    curWeights = double(curWeights);

    curWeights = log10(curWeights(curMask));

    histogram( ...
        curAx, curWeights, ...
        'BinEdges', curBinEdges, ...
        'EdgeColor', 'black');
end

curColor = cat(1, rnnConfigs.colors);
curColor = curColor(curRnnIdx, :);


xlabel(curAx, 'log10(connection strength)');
ylabel(curAx, 'Count');

Figure.config(curFig, info);
curFig.Position(3:4) = [210, 225];

% NOTE(amotta): `Figure.config` styles the background histogram. The proper
% style must thus be set here, after `Figure.config` was run.
curColors = cat(1, rnnConfigs.colors);
curColors = curColors(curRnnIds, :);
curColors = num2cell(curColors, 2);

curHists = flip(curAx.Children);
set(curHists, {'EdgeColor'}, curColors);

%% Compare RNNs in t-SNE space vs. connection strength thresholds
clear cur*;

% NOTE(amotta): Only show data points with weight percentiles below XX
curMaxPercentile = 100;

curRnnVarNames = {'threshold_percentile', 'connectivity_aa'};
curModelVarNames = {'threshPrctile', 'connectivity'};

curRnnVarNames = [curRnnVarNames, vars.rnnName];
curModelVarNames = [curModelVarNames, vars.modelName];

curMeasures = cell(numel(rnnMeasuresThresh), 1);
curTsneT = cell(numel(rnnMeasuresThresh), 1);
for curRnnIdx = 1:numel(rnnMeasuresThresh)
    curRnnType = rnnTypes(curRnnIdx);
    curRnnDataT = rnnMeasuresThresh{curRnnIdx};
    
    curRnnName = rnnTypes(1:curRnnIdx);
    curRnnName = sum(curRnnName == curRnnType);
    curRnnName = sprintf('RNN-%s-%d', curRnnType, curRnnName);
    
    % NOTE(amotta): Find checkpoint for which the binary connectome
    % (without weight threshold) comes closest to the experimentally
    % expected pairwise connectivity.
    curCkptFile = curRnnDataT(curRnnDataT.threshold_percentile == 0, :);
    curCkptFile.fracConnDiff = abs(curCkptFile.connectivity_aa - fracConn);
    curCkptFile = sortrows(curCkptFile, 'fracConnDiff', 'ascend');
    curCkptFile = curCkptFile.ckpt_file{1};
    
    curMask = strcmpi(curRnnDataT.ckpt_file, curCkptFile);
    curRnnDataT = curRnnDataT(curMask, :);
    
    % NOTE(amotta): Sort measures by increasing threshold percentile for
    % progressive visualization of data points later on.
    curRnnDataT = sortrows(curRnnDataT, 'threshold_percentile', 'ascend');
    
    % NOTE(amotta): Restirct to percentiles below XX
    curMask = curRnnDataT.threshold_percentile <= curMaxPercentile;
    curRnnDataT = curRnnDataT(curMask, :);
    curMeasures{curRnnIdx} = curRnnDataT;
    
    curRnnDataT = curRnnDataT(:, curRnnVarNames);
    curRnnDataT.Properties.VariableNames = curModelVarNames;
    
    curRnnDataT.model(:) = categorical({curRnnName});
    curRnnDataT = circshift(curRnnDataT, 1, 2);
    curTsneT{curRnnIdx} = curRnnDataT;
end

curRnnNames = cellfun( ...
    @(t) char(t.model(1)), rnnTsneT, 'UniformOutput', false);
curRnnColors = reshape({rnnConfigs.colors}, [], 1);
curRnnColors = num2cell(cell2mat(curRnnColors), 2);

curPlotConfigT = [modelPlotConfigT([], :); ...
    curRnnNames, curRnnNames, curRnnColors];

curNames = cellfun( ...
    @(data) char(data.model(1)), ...
    curTsneT, 'UniformOutput', false);
curNames = reshape(curNames, 1, []);

%% Measures
curLimX = [0, 100];
curTicksX = 0:20:100;
curTickLabelsX = arrayfun( ...
    @(x) strrep(num2str(x, '%g'), '0.', '.'), ...
    curTicksX, 'UniformOutput', false);
curTickLabelsX(2:(end - 1)) = {''};

curColors = cat(1, rnnConfigs.colors);
curFig = figure();

for curVarIdx = 1:numel(vars)
    curVar = vars(curVarIdx);
    curModelVarName = curVar.modelName;
    
    curAx = subplot(1, numel(vars), curVarIdx);
    hold(curAx, 'on');
    
    % Indicate connectivity fractions consistent with mouse barrel
    plot(curAx, repelem(fracConn, 2), curVar.yLim, 'k--');
    
    for curRnnIdx = 1:numel(curMeasures)
        curX = curMeasures{curRnnIdx}.threshold_percentile;
        curY = curMeasures{curRnnIdx}.(curVar.rnnName);
        curColor = curColors(curRnnIdx, :);
        plot(curAx, curX, curY, 'Color', curColor, 'LineWidth', 2);
    end
end

curFig.Position(3:4) = [850, 200];
curFig.Color = 'white';

curAxes = findobj(curFig, 'Type', 'Axes');
curAxes = flip(curAxes);

arrayfun(@(a) xticks(a, curTicksX), curAxes);
arrayfun(@(a) xticklabels(a, curTickLabelsX), curAxes);
arrayfun(@(a, c) set(a, 'YScale', c.yScale), curAxes(:), vars(:));
arrayfun(@(a, c) ylim(a, c.yLim), curAxes(:), vars(:));
arrayfun(@(a, c) yticks(a, c.yTicks), curAxes(:), vars(:));
arrayfun(@(a, c) ylabel(a, c.yLabel), curAxes(:), vars(:));

set(curAxes, ...
    'TickDir', 'out', ...
   {'XLim'}, {curLimX}, ...
   {'TickLength'}, {[0.05, 0.025]}, ...
   {'YMinorTick'}, {'on'});

arrayfun(@(a) axis(a, 'square'), curAxes);
xlabel(curAxes(1), 'Connection weight threshold (percentile)');

annotation( ...
    curFig, ...
    'textbox', [0, 0.9, 1, 0.1], ...
	'String', {info.filename; info.git_repos{1}.hash}, ...
    'EdgeColor', 'none', 'HorizontalAlignment', 'center');

%% Quantitative evaluation at 50th percentile
curEvalT = cat(1, curTsneT{:});
curEvalT = curEvalT(curEvalT.threshPrctile == 50, :);
curVarNames = curEvalT.Properties.VariableNames(4:end);

for curVarIdx = 1:numel(curVarNames)
    curVarName = curVarNames{curVarIdx};
    curData = curEvalT{:, curVarName};
    curData = reshape(curData, [], 2);
    curData = num2cell(curData, 1);
    
   [curMeans, curStds] = cellfun( ...
       @(d) deal(mean(d), std(d)), curData);
   [~, curPval, curBlah] = kstest2(curData{:});
   
    fprintf('%s: %.2f+-%.2f vs. %.2f+-%.2f, p=%.3f\n', ...
        curVarName, curMeans(1), curStds(1), ...
        curMeans(2), curStds(2), curPval);
end

%% tSNE
rng(6);
curTsnes = cat(1, curTsneT{:});
curTsnes = tsne(curTsnes{:, 4:end});
curTsnes = mat2cell(curTsnes, cellfun(@height, curTsneT), 2);

% Visualization using tSNE
curFig = figure();
curAx = axes(curFig);
hold(curAx, 'on');

for curDataIdx = 1:numel(curTsneT)
    curTsne = curTsnes{curDataIdx};
    curName = char(curTsneT{curDataIdx}.model(1));
    
    curColor = strcmpi(curPlotConfigT.tag, curName);
    curColor = curPlotConfigT.color(curColor, :);
    
    curSizes = repelem(12, size(curTsne, 1));
    curSizes(end) = 48;

    plot(curAx, ...
        curTsne(:, 1), curTsne(:, 2), ...
        'Color', curColor, 'LineWidth', 1.5);
    scatter(curAx, ...
        curTsne(:, 1), curTsne(:, 2), ...
        curSizes, curColor, 'filled');
end

curIsLine = arrayfun(@(p) isgraphics(p, 'Line'), curAx.Children);
[curIsLine, curSortIds] = sort(curIsLine, 'ascend');
curAx.Children = curAx.Children(curSortIds);

xlabel(curAx, 'First t-SNE dimension');
ylabel(curAx, 'Second t-SNE dimension');
axis(curAx, 'equal');

[~, curLeg] = ismember(curNames, curPlotConfigT.tag);
curLeg = curPlotConfigT.title(curLeg);

legend( ...
    flip(curAx.Children(not(curIsLine))), ...
    curLeg, 'Location', 'EastOutside');

Figure.config(curFig, info);
curFig.Position(3:4) = [950, 620];

%% Quantitative evaluation of clustering
curConnCount = unique(cellfun(@height, curTsneT));
assert(isscalar(curConnCount));

[~, ~, curTypeIds] = unique(rnnTypes);
curIntraMask = reshape(curTypeIds, [], 1) == reshape(curTypeIds, 1, []);

curThreshPrctiles = nan(1, curConnCount);
curMaxAccuracies = nan(1, curConnCount);
curDPrimes = nan(1, curConnCount);

for curConnIdx = 1:curConnCount
    curConnData = cat(1, curTsneT{:});
    curConnData = curConnData(curConnIdx:curConnCount:end, :);
    curConnThreshPrctile = unique(curConnData.threshPrctile);
    curThreshPrctiles(curConnIdx) = curConnThreshPrctile;
    curConnData = curConnData{:, 4:end};
    
    curNormConsts = prctile(curConnData, [20, 80], 1);
    curNormConsts = diff(curNormConsts, 1, 1);
    
    curDists = @(a, b) sum(abs(a - b) ./ curNormConsts, 2);
    curDists = squareform(pdist(curConnData, curDists));
    
    curUniDists = unique(curDists);
    curAccuracies = nan(size(curUniDists));
    
    % NOTE(amotta): The d' ("dee prime") sensitivity index can be
    % approximated by sqrt(2) * Z-transform(area under receiver curve).
    % See https://en.wikipedia.org/w/index.php?title=Sensitivity_index&oldid=963482225
   [~, ~, ~, curAuc] = perfcurve(curIntraMask(:), -curDists(:), true);
    curAucDPrime = sqrt(2) * norminv(curAuc);
    curDPrimes(curConnIdx) = curAucDPrime;
    
    for curDistIdx = 1:numel(curUniDists)
        curUniDist = curUniDists(curDistIdx);
        curPredMask = curDists <= curUniDist;
        
        curAccuracy = curPredMask == curIntraMask;
        curAccuracies(curDistIdx) = mean(curAccuracy(:));
    end
    
    curMaxAccuracy = max(curAccuracies(:));
    curMaxAccuracies(curConnIdx) = curMaxAccuracy;
end


% Accuracy
curFig = figure();
curAx = axes(curFig);
axis(curAx, 'square');

curLine = plot(curAx, curThreshPrctiles, curMaxAccuracies);
set(curLine, 'Color', 'black', 'LineWidth', 2);

xlim(curAx, [0, 100]); ylim(curAx, [0.5, 1]);
xlabel(curAx, 'Connection weight threshold (percentile)');
ylabel(curAx, 'Accuracy of RNN separation');

curFig.Position(3:4) = 330;
Figure.config(curFig, info);


% Sensitivity
curFig = figure();
curAx = axes(curFig);
axis(curAx, 'square');

curLine = plot(curAx, curThreshPrctiles, curDPrimes);
set(curLine, 'Color', 'black', 'LineWidth', 2);

xlim(curAx, [0, 100]); ylim(curAx, [0, 2]);
xlabel(curAx, 'Connection weight threshold (percentile)');
ylabel(curAx, 'RNN sensitivity index (d'')');

curFig.Position(3:4) = 330;
Figure.config(curFig, info);
