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

%% Configuration
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');

varNames = { ...
    'eeRelativeReciprocity', 'eiRelativeReciprocity', ...
    'ieRelativeReciprocity', 'iiRelativeReciprocity', ...
    'relativeCycles5', 'excInOutDegreeCorrelation'};
modelNames = { ...
    'er', 'exp', 'll', 'syn', ...
    'fever', 'api', 'sorn'};

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

%% Loading data
[~, modelRawT] = ...
    RNN.loadModelMeasures( ...
        modelSummaryXlsFile, modelRawXlsDir);
assert(all(ismember(modelNames, fieldnames(modelRawT))));

%% Compute normalization constants
normConsts = struct2cell(modelRawT);
normConsts = cat(1, normConsts{:});

normConsts = normConsts(:, varNames);
normConsts = table2array(normConsts);

normConsts = prctile(normConsts, [20, 80], 1);
normConsts = diff(normConsts, 1, 1);

%% Compute distances between models
expectedMinDist = nan(numel(modelNames));

for curIdxA = 1:numel(modelNames)
    curA = modelRawT.(modelNames{curIdxA});
    
    for curIdxB = 1:numel(modelNames)
        if curIdxA == curIdxB; continue; end
        
        curB = modelRawT.(modelNames{curIdxB});
        curDists = zeros(height(curA), height(curB));
        for curIdxVar = 1:numel(varNames)
            curDists = curDists + pdist2( ...
                curA.(varNames{curIdxVar}), ...
                curB.(varNames{curIdxVar}) ...
            ) / normConsts(curIdxVar);
        end
        
        expectedMinDist(curIdxA, curIdxB) = mean(min(curDists, [], 2));
    end
end

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

imagesc(log10(expectedMinDist));

for curIdxA = 1:numel(modelNames)
   [~, curIdxB] = min(expectedMinDist(curIdxA, :));
   
    plot(curAx, ...
        curIdxB + [-0.5, +0.5, +0.5, -0.5, -0.5], ...
        curIdxA + [+0.5, +0.5, -0.5, -0.5, +0.5], ...
        'black', 'LineWidth', 2);
end

curCbar = colorbar(curAx);
curCbar.Label.String = 'Expected minimum distance';
curAx.CLim = [floor(curAx.CLim(1)), ceil(curAx.CLim(2))];
curCbar.Ticks = curAx.CLim(1):curAx.CLim(2);
curCbar.TickLabels = arrayfun( ...
    @num2str, 10 .^ curCbar.Ticks, ...
    'UniformOutput', false);

curLabels = cellfun(@upper, modelNames, 'UniformOutput', false);

axis(curAx, 'ij');
axis(curAx, 'square');

xlim(curAx, [0.5, numel(modelNames) + 0.5]);
xticks(curAx, 1:numel(curLabels));
xticklabels(curAx, curLabels);

ylim(curAx, [0.5, numel(modelNames) + 0.5]);
yticks(curAx, 1:numel(curLabels));
yticklabels(curAx, curLabels);

Figure.config(curFig, info);


%% Compute ordering
[~, ordering] = sort(expectedMinDist, 2, ...
    'ascend', 'MissingPlacement', 'last');

for curIdxA = 1:numel(modelNames)
    curModelIds = ordering(curIdxA, :);
    curModelNames = strjoin(modelNames(curModelIds), ', ');
    
    fprintf( ...
        'Most similar models to %s: %s\n', ...
        modelNames{curIdxA}, curModelNames);
end
