function [STA_output]=BC_STA(loadData,filename,pixelSize,seed,Hz)
%[STA_output]=BC_STA(loadData,filename,pixelSize,seed,Hz)
%this function loads the voltage trace of the selected bipolar cell (BC)
%and perorms a response-weighted average analysis analogous to the commonn
%spike-triggered average (STA). For simplicity we refer to the analysis
%STA,even though we work with a continous voltage signal and not spikes.
%Inputs:        
%               loadData = the folder path of the white noise files
%               filename = the filename of the BC of interest
%               pixelSize = the size of a square (checker) in monitor
%               pixels (see Documentation_Data.pdf)
%               seed = seed for the random number generator
%               Hz = stimulus update rate
%
%Outpus:        STA_output = structure containing the STA, stimulus signal,
%               averaged membrane potential for calculating the output
%               nonlinearity
%
%
% Note: the function is built only for cells recorded without frozen frames.
% But additional comments are made for cells with frozen frames.
% code by HS, last modified March 2021

%% constants
%OLED dimensions for regenerating the stimulus
monitorXPixels=800; %x pixels of OLED monitor; 
monitorYPixels=600; %y pixels of OLED monitor; 
microPerPix=2.5; %one pixel=2.5 micrometer in size
STA_TimeBack=1000; %time in ms you want to go back for generating the STA (to get smaller data I use here 1000ms, in the manuscript we use 2000ms). 

%-------------------------------------------------------------------------------
%% 1.load cell
data=load([loadData,filename]); %load the file

%-------------------------------------------------------------------------------
%% 2. regenerate the stimulus
Nframes=length(data.ttls_ms); %number of frames shown (how many times the stimulus was updated) 
Nx=ceil(monitorXPixels/pixelSize); %number of squares shown X dimension
Ny=ceil(monitorYPixels/pixelSize); %number of squares shown Y dimension
nbrCheckers_total=Nx*Ny; %total number of checkers/square presented on the screen
stimulus = recreateBinaryWhiteNoiseStimulus(Nx, Ny, Nframes, seed); %the function is available under: https://github.com/gollischlab/RecreateWhiteNoiseStimuliWithMatlab 
stimulus =stimulus(:,:,1:end-1);%take the last frame out ->needed for doing the average per frames see avgPerFrame
%the output stimulus has 3 dimensions: Y,X (nbr of checkers on Y and X axis) and time. Easiest to understand the stimulus is by
%plotting for example the first stimulus frame presented on the retina:
%imagesc(stimulus(:,:,1));
%Note: the stimulus values are drawn onto the
%retina from the lower left corner, yet matlab displays it from the upper left corner.
%Therefore the y axis has to be flipped when displaying.
%set(gca,'YDir','normal');
%For me the easiest is to directly flip the stimulus y dimension (and then I dont flip the axis when displaying):
stimulus_origD = flipdim(stimulus,1);
clear stimulus
%Yet that might be personal preferences.Note: For generating the STA it does not
%matter whether it is flipped or not. However if the STA is convolved with the natural movies it matters.
%Further, if you have frozen and running frames, regenerate the stimulus separatly
%for running and frozen frames (see Documentation_Data.pdf for the corresponding seeds).

%-------------------------------------------------------------------------------
%% 3. make average per stimulus frame
membranPotAVG=avgPerFrame(data.membPot,data.ttls_ms); % weights needed for computing the STA
 
%-------------------------------------------------------------------------------
%% 4. Computing the STA
%decide for numbers of frames you want to go back in time
STA_fnbr=(ceil(STA_TimeBack/(1000/Hz))); % define the number of frame you want to go back in time
%here I prefer to work with a 2D signal of the stimulus
stimulus_2D = reshape(stimulus_origD, nbrCheckers_total, []);
tic
[~,STA_2D]=buildSTA ...
    (stimulus_2D,membranPotAVG,STA_fnbr);
toc
%if you like the STA as a 3D e.g.:
STA_3D=reshape(STA_2D,Ny,Nx, ...
    size(STA_2D,2));% Y (monitor axis), X (monitor axis) and Z (STA_fnbr)

%-------------------------------------------------------------------------------
%% 5. Make a zoom-in of the STA
%For simplification and to keep the code simple, I cut a window of around 
%100 micrometer on each side from the maximum. A zoom-in of the STA is better
%to avoid overfitting for the NL and prediction.See Methods part of the 
%manuscript and comment-section 7 below for how it is done in the manuscript.
oneCheckInMicron=pixelSize*microPerPix;
zoom_RF=ceil(100/oneCheckInMicron);

%get the maximum of the STA to make a zoom-in
[max_Color,indx_Zoom]=max(abs(STA_3D(:)));
[posY,posX,posZ]=ind2sub(size(STA_3D),indx_Zoom);
%cut a window around the STA ->make a zoom-in STA
%zoom STA -> here you take the x-y coordinates
rangeY_zoom = posY+(-zoom_RF:zoom_RF); rangeY_zoom = rangeY_zoom(rangeY_zoom>0 & rangeY_zoom<=size(STA_3D,1));
rangeX_zoom = posX+(-zoom_RF:zoom_RF); rangeX_zoom = rangeX_zoom(rangeX_zoom>0 & rangeX_zoom<=size(STA_3D,2));
STA_3D_zoom=STA_3D(rangeY_zoom,rangeX_zoom,:); %take the data from the 2000ms STA
%reshape the STA
STA2D_zoom=reshape(STA_3D_zoom,size(rangeY_zoom,2)*size(rangeX_zoom,2),size(STA_3D_zoom,3));
%get the zoom in of the stimulus signal
stimulus3D_zoom=stimulus_origD(rangeY_zoom,rangeX_zoom,:); 
%reshap to 2D
stimulus2D_zoom=reshape(stimulus3D_zoom,size(rangeY_zoom,2)*size(rangeX_zoom,2),size(stimulus3D_zoom,3));

%-------------------------------------------------------------------------------
%% 5. visualize STA: plot the frame with maximal intensity value
%I also like to plot all/multiple STA frames around the maximum.
colorLimits = [-max_Color,max_Color]; %mainly relevant if you plot multiple frames and want to keep the same colorbar
%plot
figure;
imagesc(STA_3D(:,:,posZ))
caxis(colorLimits);
colormap('gray'); %how the sitmulus was shown
axis equal; axis tight;
title(['STA best frame, filename: ' filename]) 

%-------------------------------------------------------------------------------
%% 6. add variables into the output
STA_output.STA_2D=STA_2D;
STA_output.STA2D_zoom=STA2D_zoom;
% STA_output.STA=STA; %not  normalized
STA_output.stimulus2D_zoom=stimulus2D_zoom;
STA_output.Nx=Nx; %for reshapeing to 3D structure if needed
STA_output.Ny=Ny;
STA_output.membranPotAVG=membranPotAVG;
STA_output.STA_fnbr=STA_fnbr;
STA_output.filename=filename;
STA_output.rangeX_zoom=rangeX_zoom;
STA_output.rangeY_zoom=rangeY_zoom;
STA_output.pixelSize=pixelSize;

%-------------------------------------------------------------------------------
%% 7. additional comments:
%From the zoom STA, I seperate the response into the highest-ranked spatial and temporal
%components by singular-value decomposition (see Mansucript Schreyer & Gollisch, 2021)
%example code for the SVD:
%[sp,sigma,temp] = svd(STA,'econ'); %use the non-normalized STA here!
%sp=spatial component and temp= the temporal component. 
%the highest-ranked spatial and temporal component has to be defined e.g. by manually looking
%at the components.
%From here, I fitted a 2D Gaussian to the best spatial component.
%Important for nonlinearity and prediction calculation:
%To avoid noise contributions from pixels outside the receptive field, I reduced 
%the number of elements in the STA by setting pixel values of the STA outside 
%the 3-sigma contour of the Gaussian fit to zero and re-separating the spatiotemporal 
%receptive field within this window into the highest-ranked spatial and temporal components 
%by singular-value decomposition. To compute the nonlinearity, each component 
%was normalized to unit Euclidean norm. E.g. 
%sp_bestnorm=normSTA(sp_best);
%temp_bestnorm=normSTA(temp_best);

end