% Matlab file for plotting spiking data
close all
fontname = ("Times New Roman");

num_RGC = 135;
runtime = 1500;
sample_time = 0.1;                  % in ms
sample_freq = 1000/sample_time;

RGC_param = readmatrix("Matrices/RGC_parameters.csv");
t = readmatrix("Cell responses/t.csv");
t = t(2:15001);
type = RGC_param(:, 1);

ON_RGCs = find(type == 1);
OFF_RGCs = find(type == 2);

directory = "Cell responses/RGCs_axon/";
S = dir(fullfile(directory, "*.csv"));

for k = 1:num_RGC
    F = fullfile(directory, S(k).name);
    S(k).data = csvread(F);
    S(k).data = S(k).data(4002:19001);
    [S(k).spike_amps, S(k).spike_timing] = findpeaks(S(k).data, 'MINPEAKHEIGHT', 10, 'MinPeakDistance',50);
    num_spikes(k) = size(S(k).spike_timing, 1);
end

fig = figure;
xlim([0 15001]);
% line([3000 3000], [0 70], 'Color', 'k', 'LineStyle','--');

count = 1;
for k = 1:size(ON_RGCs, 1)
    %if((num_spikes(ON_RGCs(k))) > 1)
        for i = 1:num_spikes(ON_RGCs(k))
            line([S(ON_RGCs(k)).spike_timing(i) S(ON_RGCs(k)).spike_timing(i)], [(count-0.5) (count+0.5)], 'Color', 'b');
        end
        count = count + 1;
    %end
end
ylim([0 count]);
axis off;
fig.Position = ([600 300 250 150]);

fig = figure;
xlim([0 15001]);
% line([3000 3000], [0 70], 'Color', 'k', 'LineStyle','--');

count = 1;
for k = 1:size(OFF_RGCs, 1)
    %if((num_spikes(OFF_RGCs(k))) > 1)
        for i = 1:num_spikes(OFF_RGCs(k))
            line([S(OFF_RGCs(k)).spike_timing(i) S(OFF_RGCs(k)).spike_timing(i)], [(count-0.5) (count+0.5)], 'Color', 'r');
        end
        count = count + 1;
    %end
end
ylim([0 count]);
axis off;
fig.Position = ([600 300 250 150]);

% %%
% fig = figure;
% hold on
% colormap gray
% xlim([0 1100]);
% ylim([-1 0]);
% % initial off
% x = [0 0 300 300];
% y = [-1 0 0 -1];
% c = [0 0 0 0];
% patch(x, y, c, "EdgeAlpha", 0);
% 
% 
% x = [300 300 500 500];
% y = [-1 0 0 -1];
% c = [1 1 1 1];
% patch(x, y, c, "EdgeAlpha", 0);
% 
% x = [500 500 1100 1100];
% y = [-1 0 0 -1];
% c = [0 0 0 0];
% patch(x, y, c, "EdgeAlpha", 0);
% hold off

%%
Rod = readmatrix("Cell responses/PRs/PR_lower_214.csv");
Cone = readmatrix("Cell responses/PRs/PR_lower_215.csv");
HZ = readmatrix("Cell responses/HZs/HZ_lower_29.csv");
RBC = readmatrix("Cell responses/BIPs/BIP_lower_109.csv");
ONBC = readmatrix("Cell responses/BIPs/BIP_lower_108.csv");
OFFBC = readmatrix("Cell responses/BIPs/BIP_lower_107.csv");
AC = readmatrix("Cell responses/AIIs/AII_lower_7.csv");
ONRGC = readmatrix("Cell responses/RGCs_axon/RGC_axon_005.csv");
OFFRGC = readmatrix("Cell responses/RGCs_axon/RGC_axon_079.csv");

Rod = Rod(5001:20000);
Cone = Cone(5001:20000);
HZ = HZ(5001:20000);
RBC = RBC(5001:20000);
ONBC = ONBC(5001:20000);
OFFBC = OFFBC(5001:20000);
AC = AC(5001:20000);
ONRGC = ONRGC(5001:20000);
OFFRGC = OFFRGC(5001:20000);


fig = figure;
tiledlayout(9,1, "TileSpacing", "none")

nexttile
plot(t, Rod, 'k');
xlim([0 1500]);
%ylim([-80 -40]);
axis off

nexttile
plot(t, Cone, 'm');
xlim([0 1500]);
%ylim([-80 -40]);
axis off

nexttile
plot(t, HZ, 'k');
xlim([0 1500]);
%ylim([-80 -30]);
axis off

nexttile
hold on;
line([1500 1500], [-30 -20], 'Color', 'k');
plot(t, RBC, 'k');
xlim([0 1500]);
%ylim([-40 -20]);
axis off

nexttile
hold on;
line([1500 1500], [-30 -20], 'Color', 'k');
line([1400 1500], [-20 -20], 'Color', 'k');
plot(t, ONBC, 'b');
xlim([0 1500]);
%ylim([-40 -20]);
axis off

nexttile
hold on;
line([1500 1500], [-60 -50], 'Color', 'k');
plot(t, OFFBC, 'r');
xlim([0 1500]);
ylim([-60 -30]);
axis off

nexttile
hold on;
line([1500 1500], [-60 -50], 'Color', 'k');
plot(t, AC, 'k');
xlim([0 1500]);
ylim([-70 -40]);
axis off

nexttile
hold on;
line([1500 1500], [-40 10], 'Color', 'k');
plot(t, ONRGC, 'b');
xlim([0 1500]);
ylim([-90 30]);
axis off

nexttile
hold on;
line([1500 1500], [-40 10], 'Color', 'k');
plot(t, OFFRGC, 'r');
xlim([0 1500]);
ylim([-90 30]);
axis off

fig.Position = ([400 200 200 400]);

%% 

bin_time = 10;     % ms
num_bins = runtime/bin_time;
bin_samples = bin_time/sample_time;

spike_timing_within = [];
for k = 1:num_RGC
    for i = 1:(num_bins - 1)
        count = 0;
        for j = 1:size(S(k).spike_timing, 1)
            if(S(k).spike_timing(j) > (i-1)*bin_samples && S(k).spike_timing(j) <= i*bin_samples)
                spike_timing_within(count+1) = S(k).spike_timing(j);
                count = count + 1;
            end
            total = 0;
            for l = 1:count-1
                total = total + (spike_timing_within(l+1) - spike_timing_within(l));
            end
        end
        S(k).spikecount(i) = count;
        spike_period = ((total*sample_time)/count);
        if(spike_period == 0)
            S(k).spikefreq(i) = 0;
        else
            S(k).spikefreq(i) = 1000/((total*sample_time)/count);
        end
    end
end

ON_SpikeCount = zeros((num_bins-1), 1);
for k = 1:size(ON_RGCs, 1)
    if((num_spikes(ON_RGCs(k))) > 1)
        for i = 1:(num_bins-1)
            ON_SpikeCount(i) = ON_SpikeCount(i) + S(ON_RGCs(k)).spikecount(i);
        end
    end
end

ON_SpikeRate = (ON_SpikeCount/size(ON_RGCs, 1))/(bin_time/1000);

fig = figure;
hold on;
%line([20 20], [0 70], 'Color', 'k', 'LineStyle','--');
%title("ON RGC Spike Count")
x = 1:1:(num_bins-1);
bar(x*bin_time, ON_SpikeRate,'blue', "EdgeColor", "none")
fig.Position = ([200 300 250 80]);
ylim([0, 100]);
box off

OFF_SpikeCount = zeros((num_bins-1), 1);
for k = 1:size(OFF_RGCs, 1)
    if((num_spikes(OFF_RGCs(k))) > 1)
        for i = 1:(num_bins-1)
            OFF_SpikeCount(i) = OFF_SpikeCount(i) + S(OFF_RGCs(k)).spikecount(i);
        end
    end
end

OFF_SpikeRate = (OFF_SpikeCount/size(OFF_RGCs, 1))/(bin_time/1000);

fig = figure;
hold on;
%line([20 20], [0 30], 'Color', 'k', 'LineStyle','--');
%title("OFF RGC Spike Count")
x = 1:1:(num_bins-1);
bar(x*bin_time, OFF_SpikeRate,'r', "EdgeColor", "none")
fig.Position = ([200 300 250 80]);
ylim([0, 100]);
box off

% Figures for showing reference to spontaneous activity (first 50 bins are
% spontaneous spiking rate)
ON_SpontRate = mean(ON_SpikeRate(1:50));
OFF_SpontRate = mean(OFF_SpikeRate(1:50));

% fig = figure;
% x = 1:1:(num_bins-1);
% bar(x*bin_time, ON_SpikeRate - ON_SpontRate, "FaceColor", "b", "EdgeColor", "none")
% fig.Position = ([200 300 250 80]);
% ylim([-50, 50]);
% box off
% 
% fig = figure;
% x = 1:1:(num_bins-1);
% bar(x*bin_time, OFF_SpikeRate - OFF_SpontRate, "FaceColor", "r", "EdgeColor", "none")
% fig.Position = ([200 300 250 80]);
% ylim([-50, 50]);
% box off

% Figure to show preferential activation ----------------------
Pref_act = (ON_SpikeRate-ON_SpontRate)-(OFF_SpikeRate-OFF_SpontRate);
ON_idx = (ON_SpikeRate - ON_SpontRate)/max(ON_SpikeRate);
OFF_idx = (OFF_SpikeRate - OFF_SpontRate)/max(OFF_SpikeRate);
Pref_act = ON_idx - OFF_idx;

fig = figure;
x = 1:1:(num_bins-1);
bar(x*bin_time, Pref_act/2, "FaceColor", "#7E2F8E", "EdgeColor", "none")
fig.Position = ([200 300 300 120]);
ylim([0, 1]);
yticks([0, 1]);
yline(1, "--")
box off

writematrix(ON_SpikeRate-ON_SpontRate,'ON_SpikeRate.csv') 
writematrix(OFF_SpikeRate-ON_SpontRate,'OFF_SpikeRate.csv')