%analyze gks_gna_parameter_sweep
%save the output of gks_gna_parameter_sweep (matrix_exc and matrix_res)

%load file first - takes a while

%% summary values for individual simulations
si = .025;
bool = matrix_exc>0;
st_pt = find(matrix_exc(:,end),1,'first');
end_pt = find(matrix_exc(:,end),1,'last');
mid_pt = round((end_pt - st_pt)/2)+st_pt;
max_firing = zeros(size(matrix_res));
bursting_index = max_firing;
for i = 1: size(matrix_res,2)
    for j = 1:size(matrix_res,1)
        [spike_count,ap_times_matrix] = spike_counter(matrix_res{j,i},si);
        
        max_firing(j,i) = max(spike_count);
        bool_first = ap_times_matrix < mid_pt.*si;
        bool_last = ap_times_matrix > mid_pt.*si;
        if bool_last<1
            bool_last = 1;
        end
        
        bursting_index(j,i) = sum(sum(bool_first))/sum(sum(bool_last));
        
    end
end
figure
contourf(g_na_vals,g_ks_vals,max_firing/0.6,15);
colorbar

map = [linspace(0.9,0.5,100)',linspace(.35,0.325,100)',linspace(0,0.625,100)'];

colormap(map)
caxis([5 50])
fig = gcf;
      
      % Output the figure
      filename = sprintf('param_sweep.tiff');
      print( fig, '-dtiff','-r500', filename );
figure
contourf(g_na_vals,g_ks_vals,log(bursting_index),10)
colormap(spring)


      
      
figure
plot(matrix_res{1,1}(:,15),'k-');
axis('off')
fig = gcf;
      filename = sprintf('sweep1.tiff');
      print( fig, '-dtiff','-r500', filename );


figure
plot(matrix_res{1,20}(:,15),'k-');
axis('off')
fig = gcf;
      filename = sprintf('sweep2.tiff');
      print( fig, '-dtiff','-r500', filename );
      
figure
plot(matrix_res{20,5}(:,15),'k-');
axis('off')
fig = gcf;
      filename = sprintf('sweep3.tiff');
      print( fig, '-dtiff','-r500', filename );
      
figure
plot(matrix_res{20,20}(:,15),'k-');
axis('off')
fig = gcf;
      filename = sprintf('sweep4.tiff');
      print( fig, '-dtiff','-r500', filename );
%% plot representative images



%% functions

function [spike_count,ap_times_matrix] = spike_counter(matrix_res,si)
h=size(matrix_res,2);
spike_count = zeros(1,h);
AP_threshold=0;

for i = 1:h
    k = 1;
   if (sum(matrix_res(:,i)>=AP_threshold)>0)                        % Does this response contain action potentials?
    
    boolean = matrix_res(:,i)>=AP_threshold;                        
    bob = boolean;
    for j = 2:length(boolean)
               
        if boolean(j) ~= boolean(j-1)                   % Find action potential start/end
            
            spike_count(i) = spike_count(i) + .5;       % Adds 1/2 to count for every start/end
            
        end
    end
   end
end
ap_pks_matrix = zeros(max(spike_count),size(matrix_res,2));
all_ap_locs_matrix = zeros(max(spike_count),size(matrix_res,2));
for i = 1:size(matrix_res,2)
boolean = matrix_res(:,i)>=AP_threshold;
 peak_sweep = matrix_res(:,i);
 all_aps = peak_sweep;all_aps(all_aps<AP_threshold)=AP_threshold;           %all aps in each sweep - fixed 1.4

 [all_ap_pks,all_ap_locs] = findpeaks(all_aps,'MinPeakDistance',round(1/si),'MinPeakHeight',AP_threshold);  %finds peaks and indices capped at 500hz
 
 for j = 1:length(all_ap_pks)
     ap_pks_matrix(j,i) = all_ap_pks(j);
     all_ap_locs_matrix(j,i) = all_ap_locs(j);
 end
 %ap_pks_matrix(j+1,i) = ap_pks_matrix(j,i);
  %all_ap_locs_matrix(j+1,i) = .7*Ft;    %adds final peak at end of sweep -
  %use if you want to measure ISI COV with endpoint
end
ap_times_matrix = all_ap_locs_matrix*si; %ap times in ms
% sustained_firing = (sum(ap_times_matrix>0.5,1)/0.2)';%changed to 200ms
% for i = 1:size(ap_times_matrix,2)
%     for j = 2:size(ap_times_matrix)
%         ISI_matrix(j,i) = ap_times_matrix(j,i)-ap_times_matrix(j-1,i);
%     end
% end
ap_times_matrix(ap_times_matrix==0) = NaN;
% ISI_matrix(ISI_matrix==0) = NaN;  
% ISI_matrix(ISI_matrix<0) = NaN; 
spike_count = spike_count';
end
