%% Script to create the figures for the written part of the thesis
% This script creates and modifies the figured needed for the written part
% of the thesis. The aim is to create more or less standardized figures.
% They are saved in the 'Figures' folder automatically. All figures are
% created in the order of appearance in the thesis and the numbers here
% match those in the text.
% -------------------------------------------------------------------------
% Program by: Bjarne Schultze   [last modified: 07.03.2022]
% -------------------------------------------------------------------------

%% Paragraph to create the figures for the description of the dataset
% recording example
% plot the data, catch the figure handle and the handle to the tiled layout
plotData(15,1);
fig1 = gcf;
tyldlyt = fig1.Children; 
% remove title from the tiled layout
tyldlyt.Title.String = [];
% change panel titles to A,B,C
ax = tyldlyt.Children;
ax(3).Title.String = 'A'; ax(2).Title.String = 'B'; 
ax(1).Title.String = 'C';
% improve the appearence 
enhance_figures(fig1,'KeepTitles');

% Save the figure as an image
print(fig1, "../Figures/example_recording", "-dpng", "-r1200");

% close figure window to safely reset all changes
close(fig1);

%% Paragraph to create the figures for the methods-part
% ---- for the manual analysis (examplary reactions) ----------------------
% access the data 
data13 = extract_data(13, 'no'); data13 = data13{1};
data15 = extract_data(15, 'no'); data15 = data15{1};
% set up the figure
fig2 = figure;
tiledlayout(1,3);
% plot the examplary data
nexttile    % T cell stimulated, int EPSPs 
stackedplot(data13.timeVector(108000:137000), ...
    [data13.recordingT(108000:137000,8), ...
    data13.recordingINT(108000:137000,8)], 'DisplayLabels', ...
    {["T cell","membrane potential [mV]"], ...
    ["interneuron","membrane potential [mV]"]});
xlabel("time [s]"); title("A");

nexttile    % T cell stimulated, int spikes
stackedplot(data15.timeVector(108000:137000), ...
    [data15.recordingT(108000:137000,9), ...
    data15.recordingINT(108000:137000,9)], 'DisplayLabels', ...
    ["", ""]);
xlabel("time [s]"); title("B");

nexttile    % int stimulated, T cell graded hyperpolarization
stackedplot(data13.timeVector(108000:137000), ...
    [data13.recordingT(108000:137000,9), ...
    data13.recordingINT(108000:137000,9)], 'DisplayLabels', ...
    ["", ""]);
xlabel("time [s]"); title("C");

% improve the figure appearance
enhance_figures(fig2, 'KeepTitles');

% Save the figure as an image
print(fig2, "../Figures/example_reactions", "-dpng", "-r1200");

% close figure to safely reset all figure properties
close(fig2);

% ---- for the amplitude during stimulation -------------------------------
% clipping of the test function output for explanation purposes 
% call function and catch the right figure handle
ampAtStim_funtest(32,6);
close(2)
fig3 = figure(1);
% improve the figure appearence
enhance_figures(fig3)

% get the axes handles of the tiled layout
children_tlyt = fig3.Children.Children;
ax1 = children_tlyt(3);
ax2 = children_tlyt(4);
% adjust the axes limits to get a clipping of the polts
ax1.XLim = [10, 14.5];
ax2.XLim = [10, 14.5];
% remove the legends
legend(ax1, 'off'); legend(ax2, 'off');
% label panles with A and B
ax1.Title.String = "B"; ax2.Title.String = "A";

% increase the font size of the text for both panels seperately
gchilds1 = ax1.Children;
gchilds2 = ax2.Children;
for child = 1:length(gchilds1)                  % panel 1
    if isequal(gchilds1(child).Type, 'text')
        gchilds1(child).FontSize = 12;
    end
end
for child = 1:length(gchilds2)                  % panel 2
    if isequal(gchilds2(child).Type, 'text')
        gchilds2(child).FontSize = 12;
    end
end

% Save the figure as an image
print(fig3, "../Figures/expl_current_amplitude", "-dpng", "-r1200");

% close figure window to safely reset all changes
close(fig3);

% ----- for the spike triggered average of the postsynaptic response ------
% call function and catch figure handle
spikeTrigAvg_postsyn(32);
fig4 = gcf;
% improve figure appearance 
enhance_figures(fig4);

% Save the figure as an image
print(fig4, "../Figures/spikeTrigAvg_postsyn", "-dpng", "-r1200");

% close figure window to safely reset all changes
close(fig4);

% ----- for the spike triggered average of the presynaptic potential ------
% call function and get the figure handle
spikeTrigAvg_presyn(54);
fig5 = gcf;

% adjust the figure appearance
enhance_figures(fig5)

% Save the figure as an image
print(fig5, "../Figures/spikeTrigAvg_presyn", "-dpng", "-r1200");

% close figure to safely reset all figure properties
close(fig5);

% ----- for the spike count current plot ----------------------------------
% call function and get the figure handle
plot_FI(52);
fig6 = gcf;

% adjust figure appearance
enhance_figures(fig6)

% save the figure as an image
print(fig6, "../Figures/current_spike_count", "-dpng", "-r1200");

% close figure to safely reset all figure properties
close(fig6);

% ----- for the methods valitation part -----------------------------------
% clipping of the average spike response test function
% call the function and catch the figure handle, close the others
spikeTrigAvg_postsyn_funtest(32);
close([1,2,4,5,6,7])
fig7 = gcf;

% get the axes handle
ax7 = fig7.CurrentAxes;
% modify the axes limits to get a clipping of the plot
xlim(ax7, [13.13, 13.55]);
% improve figure appearance
enhance_figures(fig7);

% Save the figure as an image
print(fig7, "../Figures/test_spikeTrigAvg_postsyn", "-dpng", "-r1200");

% close figure to safely reset all figure properties
close(fig7);

%% Create figures for the first part of the results section
% ---- for the general tendencies, most common behaviours -----------------
% get the appropriate data and set up the figure window
data1 = extract_data(14, 'no'); data1 = data1{1}; trial = 13;
fig8 = figure();
tyldlyt = tiledlayout(1,2);

selection_window = 65000:140000;

nexttile    % T cell stimulated, int EPSPs
stackedplot(data1.timeVector(selection_window), ...
    [data1.recordingT(selection_window,trial), ...
    data1.recordingINT(selection_window,trial)], 'DisplayLabels', ...
    {["T cell","membrane potential [mV]"], ...
    ["interneuron","membrane potential [mV]"]});
xlabel("time [s]"); title("A");

nexttile    % int stimulated, T cell graded hyperpolarizations
stackedplot(data1.timeVector(selection_window), ...
    [data1.recordingT(selection_window,18), ...
    data1.recordingINT(selection_window,18)], 'DisplayLabels', ...
    ["", ""]);
xlabel("time [s]"); title("B");

% improve figure appearance
enhance_figures(fig8, 'KeepTitles');

% Save the figure as an image
print(fig8, "../Figures/general_tendencies", "-dpng", "-r1200");

% close figure to safely reset all figure properties
close(fig8);

% ---- for the three main clustering criteria -----------------------------
% load the data for this criteria from the AdditionalFiles directory
load("../AdditionalFiles/decisionData.mat")
% set the font size for the histogram legends
fontSze = 12;

fig9 = figure('Position',[106,220,1325,552]);
tldlyt = tiledlayout(fig9,1,3);

% histogram for the spike amplitude
% devide the spike heights into groups resembild the groups later used for
% clusting
sh_group0 = spike_heights(spike_heights < 0);
sh_group1 = spike_heights(spike_heights < 9 & spike_heights >= 0);
sh_group2 = spike_heights(spike_heights < 30 & spike_heights >= 9);
sh_group3 = spike_heights(spike_heights >= 30);
% create the histograms with a defined bin width
nexttile
binW = 1.5;
histogram(sh_group1, 'BinLimits', [-1.5,7.5], 'BinWidth', binW, ...
    'LineWidth',1)
hold on;
histogram(sh_group2, 'BinLimits', [9,21], 'BinWidth', binW, 'LineWidth',1)
histogram(sh_group3, 'BinLimits', [33,57], 'BinWidth', binW, ...
    'LineWidth',1)
histogram(sh_group0, 'BinLimits', [-10,-8.5], 'BinWidth', binW, ...
    'LineWidth',1)
hold off;
ylim([0,28]);
% labelling
xlabel("sipke amplitude [mV]"); ylabel("count");
legend("small spikes", "medium spikes", "large spikes", "no spikes", ...
    'Location','northeast', 'FontSize',fontSze);
title("A")
% get the axes handle for the first histogram
axA = gca;

% histogram for the T cell amplitude during -2 nA int stimulation
% devide the amplitude in two groups (hyperpolarization and no
% hyperpolarization)
hypT_group0 = amp_minus2nA_stim(amp_minus2nA_stim==0 | ...
    amp_minus2nA_stim>0);
hypT_group1 = amp_minus2nA_stim(amp_minus2nA_stim < 0);
% plot the two groups in one histogram with different colors
nexttile
binW = 0.5;
histogram(hypT_group1, 'BinLimits',[-10,0], 'BinWidth',binW, ...
    'LineWidth',1, 'FaceColor', [0.3 0.6 0.85])
hold on;
histogram(hypT_group0, 'BinLimits',[0,0.5], 'BinWidth',binW, ...
    'LineWidth',1, 'FaceColor', [0.9 0.15 0.15])
hold off;
ylim([0, 28]);
% labelling
xlabel("T cell amplitude for -2 nA stimulus [mV]", ...
    'Position', [-4.697,-2.41,-1]); ylabel("count");
legend("hyperpolarization", "no hyperpolarization", ...
    'Location','northeast','FontSize',fontSze)
title("B")

% histogram for the latency of the interneuron reaction
% devide the latencies into groups resembling those used in the analysis
lat_group0 = latencies(eq(latencies,60));
lat_group1 = latencies(latencies < 5);
lat_group2 = latencies(latencies >= 5 & latencies < 60);
% plot the single groups in one histogram with different colors
nexttile
binW = 2.5;
histogram(lat_group1, 'BinLimits',[0,5], 'BinWidth',binW,'LineWidth',1, ...
    'FaceColor', [0.9 0.85 0])
hold on;
histogram(lat_group2, 'BinLimits',[5,27.5], 'BinWidth',binW, ...
    'LineWidth',1, 'FaceColor', [0 0.7 0.8])
histogram(lat_group0, 'BinLimits',[57.5,60], 'BinWidth',binW, ...
    'LineWidth',1, 'FaceColor', [0.8, 0.1 0.6])
hold off;
ylim([0,28]);
% labelling
xlabel("latency of interneuron reaction [ms]"); ylabel("count");
legend("fast reaction", "slower reaction", "no/very slow reaction", ...
    'Location','northeast','FontSize',fontSze)
title("C")
% get the axis handle for later changes
axC = gca;

% improve the appearance
enhance_figures(fig9,'KeepTitles')
fig9.Position = [106 220 1325 552];
% change the tick values for the first histogram
axA.XTick = [-10,0,10,20,30,40,50,60];
axA.XTickLabel(1) = {'no spikes'};
axA.XTickLabelRotation = -30;
axA.XLabel.Position = [22.14,-2.46,-1];
% change the tick values for the last histogram
axC.XTick = [0,10,20,30,40,50,60];
axC.XTickLabel(7) = {'no reaction'};
axC.XTickLabelRotation = -30;
axC.XLabel.Position = [30.61,-2.46,-1];

% Save the figure as an image
print(fig9, "../Figures/hist_important_features", "-dpng", "-r1200");

% close figure to safely reset all figure properties
close(fig9);

%% create figures for the clustering results within the results section
% ---- scatter plot for the bild approach ---------------------------------
% load the file containing the deanonymization information
load("../AdditionalFiles/deanonymization.mat")
% get the numbers of the cells in each package
pack_1_3 = find(eq(deanonymization.package(:,1),13));
pack_2 = find(eq(deanonymization.package(:,1),2));
pack_4_6 = find(eq(deanonymization.package(:,1),46));
pack_5 = find(eq(deanonymization.package(:,1),5));

% take the assigned cell types into account
int218 = find(eq(deanonymization.cellType(:,1),218));
int162 = find(eq(deanonymization.cellType(:,1),162));
int264 = find(eq(deanonymization.cellType(:,1),264));
int212 = find(eq(deanonymization.cellType(:,1),212));
int201203 = find(eq(deanonymization.cellType(:,1),201203));
int159 = find(eq(deanonymization.cellType(:,1),159));
int157 = find(eq(deanonymization.cellType(:,1),157));
int6061 = find(eq(deanonymization.cellType(:,1),6061));

% scatter plot with the package information included
fig10 = figure();
scatter3(spike_heights(pack_1_3),amp_minus2nA_stim(pack_1_3), ...
    latencies(pack_1_3), 'MarkerEdgeColor',[0,0,0], 'MarkerFaceColor', ...
    [0 0.66 0.85])
hold on;
scatter3(spike_heights(pack_2),amp_minus2nA_stim(pack_2),latencies(pack_2), ...
    'MarkerEdgeColor',[0,0,0], 'MarkerFaceColor',[1,0.48,0])
scatter3(spike_heights(pack_4_6),amp_minus2nA_stim(pack_4_6), ...
    latencies(pack_4_6), 'MarkerEdgeColor',[0,0,0], 'MarkerFaceColor', ...
    [0,0.87,0])
scatter3(spike_heights(pack_5),amp_minus2nA_stim(pack_5),latencies(pack_5), ...
    'MarkerEdgeColor',[0,0,0], 'MarkerFaceColor',[0.8,0,1])
hold off;
% labelling
xlabel("spike amplitude [mV]");
ylabel("T cell amplitude [mV]");
zlabel("latency [ms]");
lgnd = legend("Packages 1 & 3", "Package 2", "Packages 4 & 6", ...
    "Packages 5", 'FontSize', fontSze);

% adapt the plot style
enhance_figures(fig10);
ax = gca;
ax.CameraPosition = [228,-78,264];
lgnd.Position = [0.72,0.67,0.18,0.16];
ax.XLabel.Rotation = -5;
ax.XTick = [-10, 0, 10, 20, 30, 40, 50, 60];
ax.XTickLabel(1) = {'no spikes'};
ax.XLabel.Position = [30.4,-10.2,-5.7];
ax.XTickLabelRotation = 34;
ax.YLabel.Rotation = 34;
ax.YLabel.Position = [66.2,-4.8,-2.6];
ax.YTickLabelRotation = -6;
ax.ZTick = [0,10,20,30,40,50,60];
ax.ZTickLabel(7) = {'no reaction'};
ax.ZLabel.Position = [-14.97,-10.29,30.00];

% save the figure to a png file
print(fig10, "../Figures/scatter_blind_approach", "-dpng", "-r1200")
% close the current figure
close(fig10)

% ---- scatter plots for packages 4 to 5 ----------------------------------
% set up the figure
fig11 = figure;
tiledlayout(1,2)

% plot the data
nexttile                            % packages 4 and 6 (interneurons 60-61)
scatter3(spike_heights(int6061),amp_minus2nA_stim(int6061),...
    latencies(int6061), 'MarkerEdgeColor',[0,0,0], ...
    'MarkerFaceColor',[0.87,0.87,0.23]);
title("A")
% labelling
xlabel("spike amplitude [mV]");
ylabel("T cell amplitude [mV]");
zlabel("latency [ms]");
legend("int 60-61")

nexttile                             % package 5 (interneuron 264)
scatter3(spike_heights(int264),amp_minus2nA_stim(int264),...
    latencies(int264),'MarkerEdgeColor',[0,0,0], ...
    'MarkerFaceColor',[0,0.4,0.7]);
title("B")
% labelling
xlabel("spike amplitude [mV]");
ylabel("T cell amplitude [mV]");
zlabel("latency [ms]");
legend("int 264")

% improve the figure appearance
enhance_figures(fig11,'KeepTitles')
fig11.Position = [1,49,1536,740.8];
ax = fig11.Children.Children;

ax(2).TitleFontSizeMultiplier = 1.7;
ax(4).TitleFontSizeMultiplier = 1.7;

% align and adjust the scatter plots
ax(4).XLim = [-10,60]; ax(4).YLim = [-10,0]; ax(4).ZLim = [0,60];
ax(4).CameraPosition = [362.39,-65.15,266.87];
ax(4).XTick = [-10, 0, 10, 20, 30, 40, 50, 60];
ax(4).XTickLabel(1) = {'no spikes'};
ax(4).XLabel.Rotation = -19;
ax(4).XLabel.Position = [18.244,-7.199,-16.297];
ax(4).XTickLabelRotation = 31;
ax(4).YLabel.Rotation = 31;
ax(4).YLabel.Position = [57.94,-3.20,-6.03];
ax(4).YTickLabelRotation = -19;
ax(4).ZTick = [0,10,20,30,40,50,60];
ax(4).ZTickLabel(7) = {'no reaction'};
ax(4).ZLabel.Position = [-14.59,-10.526,30.00];
% adapt the position of the legend 
ax(3).Position = [0.395,0.741,0.07,0.03];

% align and adjust the scatter plots
ax(2).XLim = [-10,60]; ax(2).YLim = [-10,0]; ax(2).ZLim = [0,60];
ax(2).CameraPosition = [336.046,-69.564,251.027];
ax(2).XTick = [-10, 0, 10, 20, 30, 40, 50, 60];
ax(2).XTickLabel(1) = {'no spikes'};
ax(2).XLabel.Rotation = -15;
ax(2).XLabel.Position = [34.55,-9.89,-5.14];
ax(2).XTickLabelRotation = 33;
ax(2).YLabel.Rotation = 33;
ax(2).YLabel.Position = [65.95,-4.644,-0.10];
ax(2).YTickLabelRotation = -15;
ax(2).ZTick = [0,10,20,30,40,50,60];
ax(2).ZTickLabel(7) = {'no reaction'};
ax(2).ZLabel.Position = [-14.91,-10.48,30.00];
% adapt the postion of the legend
ax(1).Position = [0.89,0.763,0.062,0.03];

% Save the figure as an image
print(fig11, "../Figures/scatter_4-6", "-dpng", "-r1200");

% close figure window to safely reset all changes
close(fig11);

% ---- scatter plot for packge 2 ------------------------------------------
% plot the data
fig12 = figure;
scatter3(spike_heights(int218), amp_minus2nA_stim(int218), ...
    latencies(int218),'MarkerEdgeColor',[0,0,0], ...
    'MarkerFaceColor',[0.4,0.8,0.8])
hold on;
scatter3(spike_heights(int212), amp_minus2nA_stim(int212), ...
    latencies(int212),'MarkerEdgeColor',[0,0,0], ...
    'MarkerFaceColor',[0.8,0.2,0.1])
scatter3(spike_heights(int201203),amp_minus2nA_stim(int201203),...
    latencies(int201203), 'MarkerEdgeColor',[0,0,0],...
    'MarkerFaceColor', [0.4,0.7,0.1])
hold off;
% labelling
xlabel("spike amplitude [mV]");
ylabel("T cell amplitude [mV]");
zlabel("latency [ms]");
legend("int 218","int 212", "int 201-203")

enhance_figures(fig12)
ax = fig12.Children;

ax(2).XLim = [-10,60]; ax(2).YLim = [-10,0]; ax(2).ZLim = [0,60];
ax(2).CameraPosition = [266.6,-71.21,293.2];
ax(2).XLabel.Rotation = -8;
ax(2).XTickLabelRotation = 29;
ax(2).XTick = [-10, 0, 10, 20, 30, 40, 50, 60];
ax(2).XTickLabel(1) = {'no spikes'};
ax(2).XLabel.Position = [10.39,-4.49,-28.57];
ax(2).YLabel.Rotation = 29;
ax(2).YLabel.Position = [51.95,-0.73,-16.73];
ax(2).YTickLabelRotation = -8;
ax(2).ZTick = [0,10,20,30,40,50,60];
ax(2).ZTickLabel(7) = {'no reaction'};
ax(2).ZLabel.Position = [-14.82,-10.359,30.00];

ax(1).Position = [0.764,0.693,0.139,0.107];

% save the figure to a png file
print(fig12, "../Figures/scatter_pack2", "-dpng", "-r1200")
% close the current figure
close(fig12)

% ---- scatter plots for packages 1 and 3 ---------------------------------
% plot the data
fig13 = figure;
scatter3(spike_heights(int159), amp_minus2nA_stim(int159), ...
    latencies(int159),'MarkerEdgeColor',[0,0,0], ...
    'MarkerFaceColor',[0.1,0.8,0.2])
hold on;
scatter3(spike_heights(int157), amp_minus2nA_stim(int157), ...
    latencies(int157),'MarkerEdgeColor',[0,0,0], ...
    'MarkerFaceColor',[1,0.75,0])
scatter3(spike_heights(int162), amp_minus2nA_stim(int162), ...
    latencies(int162),'MarkerEdgeColor',[0,0,0], ...
    'MarkerFaceColor',[0.7,0.7,0.4])
hold off;
% labelling
xlabel("spike amplitude [mV]");
ylabel("T cell amplitude [mV]");
zlabel("latency [ms]");
legend("int 159","int 157", "int 162")

enhance_figures(fig13)
ax = fig13.Children;
ax(2).XLim = [-10,60]; ax(2).YLim = [-10,0]; ax(2).ZLim = [0,60];
ax(2).CameraPosition = [266.6,-71.21,293.2];
ax(2).XLabel.Rotation = -8;
ax(2).XTickLabelRotation = 29;
ax(2).XTick = [-10, 0, 10, 20, 30, 40, 50, 60];
ax(2).XTickLabel(1) = {'no spikes'};
ax(2).XLabel.Position = [10.39,-4.49,-28.57];
ax(2).YLabel.Rotation = 29;
ax(2).YLabel.Position = [51.95,-0.73,-16.73];
ax(2).YTickLabelRotation = -8;
ax(2).ZTick = [0,10,20,30,40,50,60];
ax(2).ZTickLabel(7) = {'no reaction'};
ax(2).ZLabel.Position = [-14.82,-10.359,30.00];

ax(1).Position = [0.797,0.693,0.107,0.107];

% save the figure to a png file
print(fig13, "../Figures/scatter_pack13", "-dpng", "-r1200")
% close the current figure
close(fig13)

% ---- additional figure for packages 1 and 3 -----------------------------
% get the data and set up the figure
data17 = extract_data(17,'hn'); data17 = data17{1};
fig14 = figure();
tiledlayout(1,2,'Padding','compact');

selection_window = 65000:140000;
% plot the data
nexttile            % T cell IPSPs, int stimulated
stckplt1 = stackedplot(data17.timeVector(selection_window), ...
    [data17.recordingT(selection_window,8), ...
    data17.recordingINT(selection_window,8)], 'DisplayLabels', ...
    {["T cell","membrane potential [mV]"], ...
    ["interneuron","membrane potential [mV]"]});
xlabel("time [s]"); title("A");

nexttile
scatter3(amp_1_5nA_stim(int159), amp_minus2nA_stim(int159), ...
    spike_heights(int159),'MarkerEdgeColor',[0,0,0], ...
    'MarkerFaceColor',[0.1,0.8,0.2])
hold on;
scatter3(amp_1_5nA_stim(int157), amp_minus2nA_stim(int157), ...
    spike_heights(int157),'MarkerEdgeColor',[0,0,0], ...
    'MarkerFaceColor',[1,0.75,0])
scatter3(amp_1_5nA_stim(int162), amp_minus2nA_stim(int162), ...
    spike_heights(int162),'MarkerEdgeColor',[0,0,0], ...
    'MarkerFaceColor',[0.7,0.7,0.4])
hold off;
title("B");
% labelling
xlabel("T cell amp for 1.5 nA stim [mV]");
ylabel("T cell amp for -2 nA stim [mV]");
zlabel("spike amplitude [mV]");
legend("int 159","int 157", "int 162",'Location','northeast')

% improve figure appearance
enhance_figures(fig14,'KeepTitles')
fig14.Position = [23.4,160.2,1154.4,611.2];
ax = fig14.Children.Children;

% align and adjust the scatter plot
ax(2).CameraPosition = [-5.666,-14.243,128.2412];
ax(2).XLim = [-1,0.2]; ax(2).YLim = [-2,0];
ax(2).ZTick = [-10, 0, 10, 20, 30];
ax(2).ZTickLabel = {'no spikes','0','10','20','30'};
ax(2).ZLabel.Position = [-1.103,0.133,8.155];
ax(2).XLabel.Rotation = 16;
ax(2).XLabel.Position = [-0.333,-2.147,-12.449];
ax(2).YLabel.Rotation = -33;
ax(2).YLabel.Position = [-1.119,-0.922,-11.817];
% adapt the position of the legend
ax(1).Position = [0.869,0.65,0.082,0.0985];

% Save the figure as an image
print(fig14, "../Figures/scatter_13_additional", "-dpng", "-r1200");

% close figure window to safely reset all changes
close(fig14);

%% End
% clear all variables from the workspace
clear