% Matlab file for plotting spiking data
close all
clc
clear 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(5001:20000);
    [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
        for i = 1:num_spikes(ON_RGCs(k))-1
            ON_spike_separation(k, i) = S(ON_RGCs(k)).spike_timing(i+1) - S(ON_RGCs(k)).spike_timing(i);
        end
        count = count + 1;
    %end
end

ON_spike_separation = ON_spike_separation(any(ON_spike_separation,2),:);

for k = 1:size(ON_spike_separation, 1)
    B = ON_spike_separation(k, :);
    B = B(B>0);
    ON_spikeseparation(k) = max(B);
end

ylim([0 count]);
axis off;
fig.Position = ([600 300 300 200]);

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
        for i = 1:num_spikes(OFF_RGCs(k))-1
            OFF_spike_separation(k, i) = S(OFF_RGCs(k)).spike_timing(i+1) - S(OFF_RGCs(k)).spike_timing(i);
        end
        count = count + 1;
    %end
end

OFF_spike_separation = OFF_spike_separation(any(OFF_spike_separation,2),:);

for k = 1:size(OFF_spike_separation, 1)
    B = OFF_spike_separation(k, :);
    B = B(B>0);
    OFF_spikeseparation(k) = max(B);
end

ylim([0 count]);
axis off;
fig.Position = ([600 300 300 200]);

%%
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_upper_214.csv");
Cone = readmatrix("Cell responses/PRs/PR_lower_215.csv");
HZ = readmatrix("Cell responses/HZs/HZ_lower_3.csv");
RBC = readmatrix("Cell responses/BIPs/BIP_lower_149.csv");
ONBC = readmatrix("Cell responses/BIPs/BIP_lower_179.csv");
OFFBC = readmatrix("Cell responses/BIPs/BIP_lower_181.csv");
AC = readmatrix("Cell responses/AIIs/AII_lower_5.csv");
ONRGC = readmatrix("Cell responses/RGCs_axon/RGC_axon_005.csv");
OFFRGC = readmatrix("Cell responses/RGCs_axon/RGC_axon_030.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
hold on;
line([1500 1500], [-58 -48], 'Color', 'k');
%line([1400 1500], [-40 -40], 'Color', 'k');
plot(t, Rod, 'k');
xlim([0 1500]);
%ylim([-60 -40]);
hold off;
axis off

nexttile
hold on;
line([1500 1500], [-58 -48], 'Color', 'k');
%line([1400 1500], [-58 -58], 'Color', 'k');
plot(t, Cone, 'm');
xlim([0 1500]);
%ylim([-60 -40]);
hold off;
axis off

nexttile
hold on;
line([1500 1500], [-80 -75], 'Color', 'k');
%line([1400 1500], [-60 -60], 'Color', 'k');
plot(t, HZ, 'k');
xlim([0 1500]);
%ylim([-70 -30]);
hold off;
axis off

nexttile
hold on;
line([1500 1500], [-37 -32], 'Color', 'k');
%line([1400 1500], [-53 -53], 'Color', 'k');
plot(t, RBC, 'k');
xlim([0 1500]);
%ylim([-60 -20]);
hold off;
axis off

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

nexttile
hold on;
line([1500 1500], [-50 -45], 'Color', 'k');
%line([1400 1500], [-57 -57], 'Color', 'k');
plot(t, OFFBC, 'r');
xlim([0 1500]);
%ylim([-75 -35]);
hold off;
axis off

nexttile
hold on;
line([1500 1500], [-55 -45], 'Color', 'k');
plot(t, AC, 'k');
xlim([0 1500]);
%ylim([-65 -30]);
hold off;
axis off

ONRGC_lowpass = medfilt1(ONRGC, 50);
OFFRGC_lowpass = medfilt1(OFFRGC, 50);

nexttile
hold on;
line([1500 1500], [-40 10], 'Color', 'k');
plot(t, ONRGC_lowpass, 'b');
xlim([0 1500]);
%ylim([-120 35]);
hold off;
axis off

nexttile
hold on;
line([500 500], [-40 10], 'Color', 'k');
plot(t, OFFRGC_lowpass, 'r');
xlim([0 1500]);
%ylim([-120 35]);
hold off;
axis off

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

%% 

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 300 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 300 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 ----------------------
fig = figure;
x = 1:1:(num_bins-1);
bar(x*bin_time, (ON_SpikeRate - ON_SpontRate) - (OFF_SpikeRate - OFF_SpontRate), "FaceColor", "#7E2F8E", "EdgeColor", "none")
fig.Position = ([200 300 200 120]);
ylim([0, 50]);
yticks([0, 50]);
box off

%% Calculate spike frequency

ON_spikeseparation = 10000./ON_spikeseparation;
OFF_spikeseparation = 10000./OFF_spikeseparation;

Average_ONspike = mean(ON_spikeseparation);
STD_ONspike = std(ON_spikeseparation);

Average_OFFspike = mean(OFF_spikeseparation);
STD_OFFspike = std(OFF_spikeseparation);