%Sim_Analysis.m
%Arjun Balachandar 2016
%Display and analyze data generated by running AutoSim.m, which simulates
%neuron(s) with various input parameters (see AutoSim.m). Used to visualize
%simulation data and create plots/figures

clc;
clear all;
close all;

%Load files generated after running AutoSim.m simulations
%'_ioff' refers to pre-stimulus current used to change resting Vm. If 0,
%then no pre-stimulus current used.

%Various files to load
%load('AutoSim_distim1_ioff0.mat');
%load('AutoSim_istim085_distim5_ioff0_dgA0');
%load('AutoSim_istim070_distim5_ioff0_dgA0_maxgsub15');
%load('AutoSim_distim5_ioff0_dgA0'); % ****------ CONTAINS iSTIM70 -------****
%load('AutoSim_distim5_ioff0.mat');
load('AutoSim_istim060_distim10_ioff20_dgA0_maxgsub15'); %i=60, i_off=0-20, 15x15
planeBF = @PlaneBF;

%Variable to keep constant
ind = 1; %take slice through i_stim = set-value
if num_istim > 1 %if many i_stim used by AutoSim, then will create a slice through data for each i_stim (see below)
    value = 0;
else
    value = max_istim; %if only one level of stimulus i_stim, then use this as the 'value' for creating slice (see below)
end

dx = d_gA;
dy = d_gsub;

grid_size = 15; %Width of grid (i.e. gA-gsub plot square length)

%Dimensions of grids for plots produced later
xmax = grid_size;
xmin = 0;
ymax = grid_size;
ymin = 0;

num_xTot = (max_gsub-min_gsub)/dx + 1; %for modified domain size (i.e. If grid_size is less than the maximum conductance value used (smaller plot
%than present in data)
num_x = (xmax-xmin)/dx + 1;
num_y = (ymax-ymin)/dy + 1;

%set axes
gA_axis = linspace(ymin,ymax,num_y);
gsub_axis = linspace(xmin,xmax,num_x);

%Plot Colors:
rate_c_mod = zeros((num_istim)*(num_x)*(num_y),3);
type_c_mod = zeros((num_istim)*(num_x)*(num_y),3);
rate_c = zeros((num_istim)*(num_gA)*(num_gsub),3);
type_c = zeros((num_istim)*(num_gA)*(num_gsub),3);

%For rheobase (not used here)
current_slice = zeros(num_gA,num_gsub);
rheo0_1 = zeros((num_istim)*(num_gA)*(num_gsub),3);
rheo0_2 = zeros((num_istim)*(num_gA)*(num_gsub),3);
rheo0_3 = zeros((num_istim)*(num_gA)*(num_gsub),3);
rheo0_4 = zeros((num_istim)*(num_gA)*(num_gsub),3);
rheo1_2 = zeros((num_istim)*(num_gA)*(num_gsub),3);
rheo1_3 = zeros((num_istim)*(num_gA)*(num_gsub),3);
rheo1_4 = zeros((num_istim)*(num_gA)*(num_gsub),3);
rheo2_3 = zeros((num_istim)*(num_gA)*(num_gsub),3);
rheo2_4 = zeros((num_istim)*(num_gA)*(num_gsub),3);
rheo3_4 = zeros((num_istim)*(num_gA)*(num_gsub),3);

cur_istim = -1;

cur_ind = 1;
param_array_mod = zeros((num_istim)*(num_x)*(num_y),6 + nSpikeATPs);
for i=1:num_istim*num_gA*num_gsub
    if param_array(i,2) <= grid_size && param_array(i,3) <= grid_size
        param_array_mod(cur_ind,:) = param_array(i,:);
        cur_ind = cur_ind + 1;
    end
end

%For each data-point (i.e. each neuron simulated), add pertinent data to
%comprehensive array that stores color corresponding to each neuron, based
%on firing-type or firing-rate
for i=1:(num_istim)*(num_gA)*(num_gsub)
    %Firing-Type color scheme:
    if param_array(i,4) == 0 %R = Grey
        type_c(i,1) = 0.7;
        type_c(i,2) = 0.7;
        type_c(i,3) = 0.7;
    elseif param_array(i,4) == 1 %SS = Blue
        type_c(i,1) = 0;
        type_c(i,2) = 0;
        type_c(i,3) = 1;
    elseif param_array(i,4) == 2 %DO = Green
        type_c(i,1) = 0;
        type_c(i,2) = 1;
        type_c(i,3) = 0;
    elseif param_array(i,4) == 3 %Gap = Yellow
        type_c(i,1) = 1;
        type_c(i,2) = 1;
        type_c(i,3) = 0;
    elseif param_array(i,4) == 4 %RF = Red
        type_c(i,1) = 1;
        type_c(i,2) = 0;
        type_c(i,3) = 0;
    end
    
    %Firing-rate color scheme:
    if param_array(i,5) > 0
        rate_c(i,1) = (1.0/double(param_array(i,5)));
        rate_c(i,2) = 0.4;
        rate_c(i,3) = 0.4;
    else
        rate_c(i,1) = 0.7;
        rate_c(i,2) = 0.7;
        rate_c(i,3) = 0.7;
    end
end

%Modified type_c size (if grid-size is different from gA/gsub value)
for i=1:(num_istim)*(num_x)*(num_y)
    %Firing-Type color scheme:
    if param_array_mod(i,4) == 0 %R = Grey
        type_c_mod(i,1) = 0.7;
        type_c_mod(i,2) = 0.7;
        type_c_mod(i,3) = 0.7;
    elseif param_array_mod(i,4) == 1 %SS = Blue
        type_c_mod(i,1) = 0;
        type_c_mod(i,2) = 0;
        type_c_mod(i,3) = 1;
    elseif param_array_mod(i,4) == 2 %DO = Green
        type_c_mod(i,1) = 0;
        type_c_mod(i,2) = 1;
        type_c_mod(i,3) = 0;
    elseif param_array_mod(i,4) == 3 %Gap = Yellow
        type_c_mod(i,1) = 1;
        type_c_mod(i,2) = 1;
        type_c_mod(i,3) = 0;
    elseif param_array_mod(i,4) == 4 %RF = Red
        type_c_mod(i,1) = 1;
        type_c_mod(i,2) = 0;
        type_c_mod(i,3) = 0;
    end
end

%%Keeping one variable constant, geenrate slice through data such that all
%%points in slice have that constant value of the parameter meant to be
%%constant (given by 'ind')
slize_size = num_istim*num_gA*num_gsub;
if ind == 1 %if i_stim is constant (this is used primarily)
    slice_size = num_x*num_y; %num_gsub*num_gA;
    x_ind = 2;
    y_ind = 3;
    x_size = num_gA;
    y_size = num_gsub;
    param_array_slice_contour_FR = zeros(num_y,num_x);
    param_array_slice_contour_FP = zeros(num_y,num_x);
    param_array_slice_contour_ATPave = zeros(num_y,num_x);
    param_array_slice_contour_ATP = zeros(num_y,num_x);
    param_array_slice_contour_ATP2 = zeros(num_y,num_x);
elseif ind == 2 %if gsub constant (not used)
    slice_size = num_istim*num_gA;
    x_ind = 3;
    y_ind = 1;
    x_size = num_gA;
    y_size = num_istim;
    param_array_slice_contour_FR = zeros(num_gA,num_gsub);
    param_array_slice_contour_FP = zeros(num_gA,num_gsub);
    param_array_slice_contour_ATPave = zeros(num_gA,num_gsub);
    param_array_slice_contour_ATP = zeros(num_gA,num_gsub);
    param_array_slice_contour_ATP2 = zeros(num_gA,num_gsub);
elseif ind == 3 %if gA constant (not used)
    slice_size = num_istim*num_gsub;
    x_ind = 2;
    y_ind = 1;
    x_size = num_gsub;
    y_size = num_istim;
    param_array_slice_contour_FR = zeros(num_gA,num_gsub);
    param_array_slice_contour_FP = zeros(num_gA,num_gsub);
    param_array_slice_contour_ATPave = zeros(num_gA,num_gsub);
    param_array_slice_contour_ATP = zeros(num_gA,num_gsub);
    param_array_slice_contour_ATP2 = zeros(num_gA,num_gsub);
end

%note: slice_size is modified above if modified domain size used
param_array_slice = zeros(slice_size, 6 + nSpikeATPs);
temp_slice = zeros(num_y, num_x); 
prev_slice = zeros(num_y, num_x);
rheo_slice = [];
rheo_colour = [];
color_array_slice = zeros(slice_size, 3);

cur_ind = 1;
cur_ind2 = 1;
cur_istim = min_istim;
prev_fp = 0;
cur_fp = 0;


for i=1:(num_istim)*(num_gA)*(num_gsub)

    %Calculation of rheobase (not used here)
%     if param_array(i,1) == cur_istim
%         temp_slice( (param_array(i,3) - min_gA)/d_gA + 1, (param_array(i,2) - min_gsub)/d_gsub + 1) = param_array(i,4);
%     else
%         for j=0:num_gA-1
%             for k=0:num_gsub-1
%                 if temp_slice(j+1,k+1) > prev_slice(j+1,k+1)
%                     prev_fp = prev_slice(j+1,k+1);
%                     cur_fp = temp_slice(j+1,k+1);
%                     rheo_slice(cur_ind2,1) = min_gA + d_gA*j;
%                     rheo_slice(cur_ind2,2) = min_gsub + d_gsub*k;
%                     rheo_slice(cur_ind2,3) = cur_istim; %current i_stim
%                     rheo_slice(cur_ind2,4) = temp_slice(j+1,k+1);
%                     
%                     if temp_slice(j+1,k+1) == 0 %R = Grey
%                         rheo_colour(cur_ind2,1) = 0.7;
%                         rheo_colour(cur_ind2,2) = 0.7;
%                         rheo_colour(cur_ind2,3) = 0.7;
%                     elseif temp_slice(j+1,k+1) == 1 %SS = Blue
%                         rheo_colour(cur_ind2,1) = 0;
%                         rheo_colour(cur_ind2,2) = 0;
%                         rheo_colour(cur_ind2,3) = 1;
%                     elseif temp_slice(j+1,k+1) == 2 %DO = Green
%                         rheo_colour(cur_ind2,1) = 0;
%                         rheo_colour(cur_ind2,2) = 1;
%                         rheo_colour(cur_ind2,3) = 0;
%                     elseif temp_slice(j+1,k+1) == 3 %Gap = Yellow
%                         rheo_colour(cur_ind2,1) = 1;
%                         rheo_colour(cur_ind2,2) = 1;
%                         rheo_colour(cur_ind2,3) = 0;
%                     elseif temp_slice(j+1,k+1) == 4 %RF = Red
%                         rheo_colour(cur_ind2,1) = 1;
%                         rheo_colour(cur_ind2,2) = 0;
%                         rheo_colour(cur_ind2,3) = 0;
%                     end
%                     cur_ind2 = cur_ind2 + 1;
%                 end
%             end
%         end
%         prev_slice = temp_slice;
%         temp_slice = zeros(num_gA, num_gsub);
%         cur_istim = cur_istim + d_istim;
%     end
    
    %***Create 'slice' (i.e. all data such that parameter-of-interest
    %equals***
    %'value' (inputted by user) -> e.g. all points where i_stim = value)
    if param_array(i,ind) == value && param_array(i,2) <= grid_size && param_array(i,3) <= grid_size
        for j=1:(6+nSpikeATPs)
            param_array_slice(cur_ind,j) = param_array(i,j);
        end
        color_array_slice(cur_ind,1) = type_c(i,1);
        color_array_slice(cur_ind,2) = type_c(i,2);
        color_array_slice(cur_ind,3) = type_c(i,3);
        cur_ind = cur_ind + 1;
    end
end

gridX = (1:x_size);
gridY = (1:y_size);

rheo_01X = [];
rheo_01Y = [];
rheo_01 = [];

cur_ind = 1;
for i=1:size(rheo_slice)
    if rheo_slice(i,4) == 1
        rheo_01X(cur_ind) = rheo_slice(i,1);
        rheo_01Y(cur_ind) = rheo_slice(i,2);
        rheo_01(cur_ind,cur_ind) = rheo_slice(i,3);
        cur_ind = cur_ind + 1;
    end
end

%%Display figures:

%3D-scatter plot of all points/neurons (i.e. i_stim vs. gA vs. gsub)
if num_istim > 1 
    figure('name','3D Scatter Plot - Firing Pattern');
    scat3 = scatter3(param_array_mod(:,2),param_array_mod(:,3),param_array_mod(:,1),10,type_c_mod,'filled');
    ylabel('gA');
    xlabel('gsub');
    zlabel('istim');
end

%2D-scatter plot (slice) of all points/neurons with i_stim = value (set by user) -> (i.e. i_stim vs. gA vs. gsub)
figure('name','Scatter Plot Slice - Firing Pattern');
scat = scatter(param_array_slice(:,x_ind),param_array_slice(:,y_ind),10,color_array_slice,'filled');
title(['istim = ',int2str(value),', ioff = ',int2str(i_off)]);
xlabel('gsub');
ylabel('gA');

%Same as 2D scatter plot above, but using contourf (commented here since scatter plot conveys same information):
% figure('name','Contour Plot Slice - Firing Pattern');
% [aC aC] = contourf(param_array_slice_contour_FP,5);
% title(['istim = ',int2str(value),', ioff = ',int2str(i_off)]);
% xlabel('gsub');
% ylabel('gA');
% set(aC,'LineStyle','none');
% set(gca, 'XTick', 1:x_size); % Change x-axis ticks
% set(gca, 'XTickLabel', gsub_axis); % Change x-axis ticks labels.
% set(gca, 'YTick', 1:y_size); % Change x-axis ticks
% set(gca, 'YTickLabel', gA_axis); % Change x-axis ticks labels.

%Plot firing-rate slice using contourf
% figure('name','Contour Plot Slice - Firing Rate');
% [hC hC] = contourf(param_array_slice_contour_FR,255);
% set(hC,'LineStyle','none');
% title(['istim = ',int2str(value),', ioff = ',int2str(i_off)]);
% xlabel('gsub');
% ylabel('gA');
% set(gca, 'XTick', 1:x_size); % Change x-axis ticks
% set(gca, 'XTickLabel', gsub_axis); % Change x-axis ticks labels.
% set(gca, 'YTick', 1:y_size); % Change x-axis ticks
% set(gca, 'YTickLabel', gA_axis); % Change x-axis ticks labels.
% colorbar;


%plot rheo-base (not used here)
if size(rheo_slice) > 0
    figure('name','3D Scatter Plot - Rheobase');
    scat_rheo3D = scatter3(rheo_slice(:,1),rheo_slice(:,2),rheo_slice(:,3),20,rheo_colour,'filled');
    xlabel('gA');
    ylabel('gsub');
    zlabel('istim');
end