function Analyze_Spatial_Contrast_Adaptation_Stimulus(varargin)
%
%
%%% Analyze_Spatial_Contrast_Adaptation_Stimulus %%%
%
% This function analyzes the response of the ganglion cells to the spatial
% contrast adaptation stimulus. This is a supplement code for the Khani & Gollisch
% (2017) study about the diversity of spatial contrast adaptation in the mouse retina.
% This function accompanys dataset of Khani and Gollisch (2017). The dataset is available at:
% https://gin.g-node.org/gollischlab/Khani_and_Gollisch_2017_RGC_spiketrains_for_spatial_contrast_adaptation
% This function first loads the data in a format that is provided in by the
% supplementary database of the study. It then generate PSTHs, spike-triggered
% average (STA) and nonlinearity functions for the selected cell.
%
% ===============================Inputs====================================
%
%   datapath        :   path to the folder in the database.
%   stimulusnumber  :   number to select correct stimulus.
%   cellnumber      :   an example cell number from the list of good cells.
%   clusternumber   :   a cluster number for the selected cell.
%
%================================Output====================================
%
%   PSTH, STAs and nonlinearities for the selected cell.
%

%--------------------------------------------------------------------------------------------------%
%----------                              main-function                                   ----------%
%--------------------------------------------------------------------------------------------------%

% first load the data for an example cell
[dp, stimulusnumber, cellnumber, clusternumber] = select_example_cell(varargin{:});
% analyze the data
ca_data     =       sca_anylsis(dp, stimulusnumber, cellnumber, clusternumber);
% plotting
plot_sca(ca_data, cellnumber, clusternumber);

end

%--------------------------------------------------------------------------------------------------%
%----------                               sub-functions                                  ----------%
%--------------------------------------------------------------------------------------------------%

function expdata = sca_anylsis(dp, stimulusnum, cellnum, clusternum,varargin)
%
expdata             =   load_data_for_selected_cell(dp, stimulusnum, cellnum, clusternum);
% set all the required parameters
para                =   expdata.stimpara;
para.refreshrate    =   60;
para.upsampleratio  =   ceil(1e3/para.refreshrate);
para.numNLbins      =   40;
para.meanIntensity  =   0.5;
para.statime = fliplr(linspace(0,(1e3/60*para.Nblinks)*(para.numNLbins/para.upsampleratio),para.numNLbins));
Nruns = floor((length(expdata.frametimes)-1)/(para.Nblinks*para.switchFrames));
para.Nframes = Nruns*para.Nblinks*para.switchFrames;
para.Nrepeats           =       floor(Nruns/2);   % number of runs per trial;
para.binlength          =       250 / 1e3;        % 250 ms each bin
% for the PSTH
para.low_contduration   =    para.switchFrames ./ para.refreshrate;
para.high_contduration  =    para.switchFrames ./ para.refreshrate;
para.psthtime_low   = linspace(0, para.low_contduration ,(para.low_contduration / para.binlength)+1);
para.psthtime_high  = linspace(0, para.high_contduration ,(para.high_contduration / para.binlength)+1);
para.trialduration  = para.low_contduration + para.high_contduration;
para.psthtimelh     = linspace(0, para.trialduration,(para.trialduration / para.binlength));

% loading clusters and frametimings
ft              =       expdata.frametimes;    % work in seconds
ftimes          =       ft(1:para.Nframes);
% get frame times per trial
ftimes          =       reshape(ftimes,para.switchFrames,Nruns)';
ftimes_high     =       ftimes(2:2:end,:);      ftimes_high = ftimes_high(1:para.Nrepeats,:);
ftimes_low      =       ftimes(1:2:end,:);      ftimes_low  = ftimes_low(1:para.Nrepeats,:);

% getting the stimulus right!
stimulus = Spatial_Contrast_Adaptation_Stimulus(para.seed,para.Nframes,para.switchFrames,Nruns, para.meanIntensity);

% binning spikes
spkbins         =       (histc(expdata.spiketimes,ft(1:para.Nframes))); %#ok
spkbins         =       reshape(spkbins,para.switchFrames,Nruns)';
spkbins_high    =       spkbins(2:2:end,:);      spkbins_high = spkbins_high(1:para.Nrepeats,:);
spkbins_low     =       spkbins(1:2:end,:);      spkbins_low  = spkbins_low(1:para.Nrepeats,:);

% setting output
expdata.ftimes.high         =       ftimes_high;
expdata.ftimes.low          =       ftimes_low;
expdata.stimulus            =       stimulus;
expdata.stimpara            =       para;

% calculating rasters and psth
expdata     =       calc_rasters_psth(expdata);

% calculating STAs and nonlinearities for each contrast and location
[locX_high.sta,locX_high.nlx,locX_high.nly,locX_high.time]  = calculate_sta_nonlinearity(spkbins_high, expdata.stimulus.locX_high, para);
[locX_low.sta,locX_low.nlx,locX_low.nly,locX_low.time]      = calculate_sta_nonlinearity(spkbins_low, expdata.stimulus.locX_low, para);
[locY_high.sta,locY_high.nlx,locY_high.nly,locY_high.time]  = calculate_sta_nonlinearity(spkbins_high, expdata.stimulus.locY_high, para);
[locY_low.sta,locY_low.nlx,locY_low.nly,locY_low.time]      = calculate_sta_nonlinearity(spkbins_low, expdata.stimulus.locY_low, para);

% setting output
expdata.locXhigh            =       locX_high;
expdata.locX_low            =       locX_low;
expdata.locY_high           =       locY_high;
expdata.locY_low            =       locY_low;

end

%--------------------------------------------------------------------------------------------------%

function expdata = calc_rasters_psth(expdata)

fthigh          =       expdata.ftimes.high;
ftlow           =       expdata.ftimes.low;
spks            =       expdata.spiketimes;
para            =       expdata.stimpara;
[ras_high, ras_low]     =   deal(cell(para.Nrepeats,1));
delay           =       50/1e3; % ignore first spikes after contrast switch

for ii = 1:para.Nrepeats
    
    ras_high{ii}    = spks(and(spks >= fthigh(ii,1),spks <= fthigh(ii,end))) - fthigh(ii,1);
    ras_low{ii}     = spks(and(spks >= (ftlow(ii,1)+delay),spks <= ftlow(ii,end))) - (ftlow(ii,1)+delay);
    
end
rasters_high    =       CelltoMatUE(ras_high);
rasters_low     =       CelltoMatUE(ras_low);

psth_high       =   mean(histc(rasters_high',para.psthtime_high),2)*(1/para.binlength);   %#ok
psth_high       =   psth_high(1:end-1);     psth_high(end) = psth_high(end-1);
para.psthtime_high  =   para.psthtime_high(1:end-1);
psth_low        =   mean(histc(rasters_low',para.psthtime_low),2)*(1/para.binlength);   %#ok
psth_low        =   psth_low(1:end-1);      psth_low(end) = psth_low(end-1);
para.psthtime_low   =   para.psthtime_low(1:end-1);

expdata.rasters_high        =       rasters_high;
expdata.rasters_low         =       rasters_low;
expdata.psth_high           =       psth_high';
expdata.psth_low            =       psth_low';
expdata.psthlh              =       [psth_low; psth_high]';

end

function [staNorm,nlx,nly,timeVec] = calculate_sta_nonlinearity(runningbin, stimulus, para)

toHzratio       =       para.Nblinks/para.refreshrate; % normalize the response to Hz

rnstimEn    =   zeros(para.Nrepeats,para.numNLbins, para.switchFrames-para.numNLbins+1);
for itrial = 1 : para.Nrepeats
    stimtrial        =     stimulus(itrial,:);
    hankelMat        =     hankel(1:para.numNLbins,para.numNLbins:para.switchFrames);
    rnstimEn(itrial,:,:)    =   stimtrial(hankelMat);
end
runningbin  =   reshape(permute(runningbin(:,para.numNLbins:end,:),[1 3 2]),1,[]);
rnstimEn    =   reshape(permute(rnstimEn, [2 1 3]), para.numNLbins,[]);

sta         =   bsxfun(@rdivide,runningbin*rnstimEn',sum(runningbin,2));
staNorm     =   bsxfun(@times,sta,1./sqrt(sum(sta.^2,2)));


timeVec     =   (-(para.numNLbins-1):1:0)*toHzratio; %in seconds

runningGenerators   =   staNorm * rnstimEn;

[nly,nlx]   =   calculate_nonlinearity(runningGenerators', runningbin', para.numNLbins, toHzratio);


end

%--------------------------------------------------------------------------------------------------%


function [ values, centers,sevalues] = calculate_nonlinearity(linearOutput, spikes, nbins, dt)

qtStep          = 1/nbins;
qtEdges         = 0:qtStep:1;
genEdges        = quantile(linearOutput, qtEdges);

if all(isnan(genEdges))   % this is for cases without spikes where everything is NaNs
    values = zeros(size(genEdges));
    sevalues = zeros(size(genEdges));
else
    [a, spikebins]  = histc(linearOutput,genEdges); %#ok
    
    values          = accumarray(spikebins, spikes', [nbins+1 1], @sum);
    sdvalues        = accumarray(spikebins, spikes', [nbins+1 1], @std);
    
    values          = values./a/dt;
    sevalues        = sdvalues./sqrt(a)/dt;
    
end
values(end)     = [];
sevalues(end)   = [];
centers         = genEdges(1:nbins)+diff(genEdges)/2;
values          = transpose(values);
sevalues        = transpose(sevalues);

end

%--------------------------------------------------------------------------------------------------%

function [dp, stimulusnumber, cellnumber, clusternumber] = select_example_cell(varargin)
% first parse the user inputs
p = inputParser();      % check the user options.
p.addParameter('datapath', [],  @(x) ischar(x) || isstring(x));
p.addParameter('stimulusnumber', [], @isnumeric);
p.addParameter('cellnumber', [], @isnumeric);
p.addParameter('clusternumber', [], @isnumeric);
p.addParameter('help', false,  @(x) islogical(x) || (isnumeric(x) && ismember(x,[0,1])));

p.parse(varargin{:});
% defualt parameters
defpara     =    p.Results;

if isempty(defpara.datapath)
    dp                  =       uigetdir('Select an experiment folder:');
else
    dp                  =       defpara.datapath;
end

if isempty(defpara.cellnumber) || isempty(defpara.clusternumber)
    listofgoodcells     =       load([dp,filesep,'list_of_good_cells.txt'],'-ascii');
    cellabels   =   sprintfc('cell %g, cluster %g',listofgoodcells);
    if numel(cellabels) > 1
        idx         =   listdlg('PromptString',{'Select an example ganglion cell:'},...
            'SelectionMode','single','ListString',cellabels,'listsize',[250 250]);
    else
        idx         =   1;
    end
    cellnumber          =       listofgoodcells(idx,1);
    clusternumber       =       listofgoodcells(idx,2);
else
    cellnumber      =   defpara.cellnumber;
    clusternumber   =   defpara.clusternumber;
end

if isempty(defpara.stimulusnumber)
    stimnames           =       importdata([dp,filesep,'stimuli_names.txt']);
    canames     =   stimnames(~(cellfun('isempty',(strfind(stimnames,'Contrast_Adaptation')))));
    
    if numel(canames) > 1
        idx     = listdlg('PromptString',{'Select an example dataset:'},...
            'SelectionMode','single','ListString',canames,'listsize',[450 100]);
        canames     =   canames{idx};
    end
    stimulusnumber      =       str2double(extractBefore(canames,'_'));
else
    stimulusnumber  =       defpara.stimulusnumber;
end

end

%--------------------------------------------------------------------------------------------------%

function expdata = load_data_for_selected_cell(dp, stimnum, cellnum, clusternum)
%
ftfolder        =       [dp,filesep,'frametimes'];
if ~exist(ftfolder,'dir'),      ftfolder    =   dp;     end
rasfolder       =       [dp,filesep,'spiketimes'];
if ~exist(rasfolder,'dir'),      rasfolder    =   dp;     end
paramfolder     =       [dp,filesep,'stimulusparameters'];
if ~exist(paramfolder,'dir'),      paramfolder    =   dp;     end

ftnames         =       dir([ftfolder, filesep, num2str(stimnum,'%02g'),'*_frametimings.txt']);
rastername      =       [rasfolder, filesep, sprintf('%d_SP_C%d%02d.txt', stimnum, cellnum, clusternum)];

warning('off','all');   % this is supress the warning about long names, hopefully nothing happen in between these lines!
expdata.spiketimes      =       load(rastername, '-ascii');
expdata.frametimes      =       load([ftfolder,filesep,ftnames.name], '-ascii');

stimpath                =       dir([paramfolder,filesep,num2str(stimnum,'%02g'),'*_parameters.txt']);

% now reading the txt file with parameters
fid             =       fopen([paramfolder,filesep,stimpath.name]);
tline           =       fgetl(fid);
while ischar(tline)
    tline       =       fgetl(fid);
    if tline    ==      -1,     break;      end
    fn          =       extractBefore(tline,' = ');
    if isempty(fn),     continue;           end
    val         =       extractAfter(tline,' = ');
    % to convert to double
    if ~isnan(str2double(val))
        val     =       str2double(val);
    end
    % to convert to logical
    if strcmp(val,'true'),        val = true;           end
    if strcmp(val,'false'),       val = false;          end
    % store the values in a structure
    stimpara.(fn)       =       val;
end
fclose(fid);

expdata.stimpara    =       stimpara;
expdata.name        =       extractBefore(ftnames.name,'_frametimings.txt');
expdata.paraname    =       extractBefore(stimpath.name,'_parameters.txt');
if strcmp(expdata.name , expdata.paraname)
    expdata         =       rmfield(expdata,'paraname');
end
warning('on','all');

end

%--------------------------------------------------------------------------------------------------%

function OutputMat = CelltoMatUE (input_cell)
%
inptsize        =       size(input_cell);
maxLength       =       max(cellfun('length',input_cell));
if inptsize(2)  >       inptsize(1)
    OutputMat   = cell2mat(cellfun(@(x)cat(1,x,nan(maxLength-length(x),1)),input_cell,'un',0));
else
    OutputMat   = cell2mat(cellfun(@(x)cat(2,x',nan(1,maxLength-length(x))),input_cell,'un',0));
end

end

%--------------------------------------------------------------------------------------------------%
%----------                             plotting modules                                 ----------%
%--------------------------------------------------------------------------------------------------%

function plot_sca(ca_data, cellnum, clusternum)

% putting data together to make plotting easier
para    =   ca_data.stimpara;
sta     =   [ca_data.locXhigh.sta;ca_data.locX_low.sta;ca_data.locY_high.sta;ca_data.locY_low.sta]';
statime = - [ca_data.locXhigh.time;ca_data.locX_low.time;ca_data.locY_high.time;ca_data.locY_low.time]'*1e3;
nlx     =   [ca_data.locXhigh.nlx;ca_data.locX_low.nlx;ca_data.locY_high.nlx;ca_data.locY_low.nlx]';
nly     =   [ca_data.locXhigh.nly;ca_data.locX_low.nly;ca_data.locY_high.nly;ca_data.locY_low.nly]';

figure('pos',[600 50 900 950],'color','w');
% get plots title and legend labels
titr    =   {'Location X', 'Location Y'};
if para.twofixed
    l   =   {'high-low','low-low'};
    s   =   'Spatial contrast adaptation';
else
    l   =   {'high-low','low-high'};
    s   =   'Spatial contrast adaptation (alternating stimulus)';
end

% first plotting psths
subplot(3,2,[1,2])
bar(para.low_contduration+para.psthtime_high(1:end-1),ca_data.psth_high);
hold on;
bar(para.psthtime_low(1:end-1),ca_data.psth_low);
pbaspect([3 1 1]);          box off;
xlabel('Time (sec)');       ylabel('Rate (Hz)');
xlim([0 para.low_contduration+para.high_contduration]);
xticks([0 para.low_contduration para.low_contduration+para.high_contduration]);
yAx = get(gca,'YLim');      yAx = ceil(yAx(2)/2)*2;       yticks([0 yAx/2 yAx]);
legend(l,'Location','northwest');               legend boxoff;

% plotting STAs and nonlinearities
for ii = 1:2
    
    subplot(3,2,ii+2)
    plot(statime(:,2*(ii-1)+1),sta(:,2*(ii-1)+1),statime(:,2*ii),sta(:,2*ii),'LineWidth',2);
    pbaspect([1.2 1 1]);    box off;
    axis([0 600 -0.5 0.5]);
    title(['STA ',titr{ii}]);
    xlabel('Delay (sec)');   ylabel('Filter (a.u.)');
    if ii==2, legend(l,'Location','southeast');      legend boxoff; end
    
    subplot(3,2,ii+4)
    plot(nlx(:,2*(ii-1)+1),nly(:,2*(ii-1)+1),nlx(:,2*ii),nly(:,2*ii),'LineWidth',2);
    pbaspect([1.2 1 1]);     box off;    ax = gca;
    yAx = ceil(ax.YLim(2)/2)*2;
    axis([-3.1 3.1 0 yAx]);
    title(['Nonlinearity ',titr{ii}]);
    xlabel('Input');   ylabel('\Delta Rate (Hz)');
    %ax.XAxisLocation = 'origin';    ax.YAxisLocation = 'origin';
    yticks(0:yAx/2:yAx);
    xticks(-3:1:3);
    
end
% mega title on top
annotation('textbox', [0.1, 0.95, 0.8, 0], 'string', [s, ' properties of cell ',...
    num2str(cellnum),', cluster ',num2str(clusternum)],'FontSize',16,'EdgeColor','none',...
    'HorizontalAlignment','center','VerticalAlignment','middle');

end
%--------------------------------------------------------------------------------------------------%