function [naturalMovie_output]=BC_predNaturalMovies(loadData,filename,movieSignal,STA_output,NL_output,Hz)
%[naturalMovie_output]=BC_predNaturalMovies(loadData,filename,movieSignal,STA_output,NL_output,Hz)
%this function shows the prinicpal calculation for the prediction of the 
%response to the natural movie, from the linear-nonlinear model (LN-model)
%computed with white noise stimulus.
%Inputs:
%               loadData = the folder path of the natural movie files
%               filename = the filename of the BC of interest
%               movieSignal = the preprocessed natural movie signal the
%               output of the BC_getNaturalMovies.m. Important that the
%               movie signal has the same dimensions as the STA_output
%               STA_output =  output of the BC_STA function; containing the
%               STA (filter) of the white noise stimulus  
%               NL_output =  output of the BC_NL function; containing the
%               output nonlinearity of the white noise stimulus  
%               Hz = movie update rate given in excel file
%
%Outpus:        naturalMovie_output = structure containing the output of the
%               prediction
%
% code by HS last modiefied March 2021

%% Constants
findBreak=500; %the movie is repeated multiple times. Between the movies there is 
%a break/gray screen of >0.5s. To find the start of each movie trial, you find the gray screen between the movies

%-------------------------------------------------------------------------------
%% 1. cut the natural movie signal to the size of the filter
%in the BC_STA.m function, the STA was cut. Here the natural movie signal
%has to be cut to the same dimensions.
movieSignal_zoom=movieSignal(STA_output.rangeY_zoom,STA_output.rangeX_zoom,:); 
%make the signal 2D for prediction
movieSignal2D_zoom=reshape(movieSignal_zoom,size(movieSignal_zoom,1)*...
    size(movieSignal_zoom,2),size(movieSignal_zoom,3));

%-------------------------------------------------------------------------------
%% 2. compute the prediction
%% 2.1 convolution
%normalize the filter
linearFiler=normSTA(STA_output.STA2D_zoom);
%convolve the natural movie stimulus sequence  with the
%normalized STA.
stim2pred=filter2(linearFiler,movieSignal2D_zoom,'valid');
%%IMPORTANT: if you have a best spatial and temporal component as described in Point 7
%in the function BC_STA.m (see Method of Manuscript), this part hast to be done in two lines, e.g.:
%NL_PredXaxis2=filter2(sp_bestnorm,movieSignal2D_zoom,'valid'); %first
%convolve the spatial dimension
%generator_signal=filter2(temp_bestnorm,NL_PredXaxis2,'valid'); %then the
%temporal dimension
%%Also note, to avoid noise contributions from pixels outside the receptive
%%field,it is important to cut the STA into a smaller region around the maximum pixel, this
%%gives a much cleaner NL-> see BC_STA.m and Methods in Schreyer&Gollisch (2021).

%% 2.2 compute the prediction through linear interpolation
NL_x=NL_output.bins.NL_xbin;
NL_y=NL_output.bins.NL_ybin;
%example of how to make the prediction through linear interpolation
predMemP_y=nan(1,size(stim2pred,2));
for kk=1:size(stim2pred,2)
    nbrTemp=stim2pred(kk); %first x-axis value to predict
    [xIndxS,xIndxL]= getIndex(NL_x,nbrTemp,1); %get those values that are on left and rigth
    xValueS=NL_x(xIndxS);
    xValueL=NL_x(xIndxL);
    %do the linear interpolation
    [polyCoeff] = polyfit([xValueS;xValueL],[NL_y(xIndxS);...
        NL_y(xIndxL)],1);
    predMemP_y(kk)=polyCoeff(1)*nbrTemp+polyCoeff(2);
end

%-------------------------------------------------------------------------------
%% 4. compute the average response to the move
%% 4.1 load movie response
data=load([loadData,filename]); %load the file

%sort the different pulses for the different trials
diffttls=diff(data.ttls_ms);
nbrTrials=find(diffttls>findBreak); %to find out which pulse correspond to a new trial (gray break)
nbrFrames=diff(nbrTrials);
if isempty(nbrFrames)
    nbrFrames=nbrTrials;
end
%I get the membrane potential per trial & frames per trial
reOrderFrames=nan(nbrFrames(1),length(nbrTrials));
membranePerTrial_origRes=cell(1,length(nbrTrials));
for kk=1:length(nbrTrials)
    if kk==1 %first trial
        reOrderFrames(:,kk)=data.ttls_ms(1:nbrTrials(kk));
        membranePerTrial_origRes{kk}=data.membPot(reOrderFrames(1,kk):reOrderFrames(end,kk)-1); %-1 because the last frame is the start of the gray screen
    else
        reOrderFrames(:,kk)=data.ttls_ms(nbrTrials(kk-1)+1:nbrTrials(kk)); %end of previous trial +1 (starts at 999)
        membranePerTrial_origRes{kk}=data.membPot(reOrderFrames(1,kk):reOrderFrames(end,kk)-1);
    end
end
%I only include trials here that here fully finished until the break!

%------------------------------------------------------------------------------------------
%% 4.2 Natural Movie Build average per frame
%make one value for the membrane potential between two stimulus signal (two
%frames), where the same image was shown
memPperTrial_Movie=avgPerFrameNatMovies(data.membPot,...
    reOrderFrames,nbrFrames(1),Hz);

memP_AVGMovie=mean(memPperTrial_Movie,2); 

%------------------------------------------------------------------------------------------
%% 5 compute the pearson correlation coefficient
%pearson correlation
corrcoeff_temp=corrcoef(memP_AVGMovie(STA_output.STA_fnbr:end),predMemP_y); 
pearCorrSqrt=corrcoeff_temp(2)^2;
%not with this method the prediction value is a bit lower than what we report in the
%manuscript. To get the same values as in the mansucript, you have to
%cut 3sigma around the RF and perform the SVD (see also Manuscript and
%BC_STA, BC_NL)

%-------------------------------------------------------------------------------
%% 6. plot the output
figure;
%plot the first frame of the move
subplot(2,1,1)
imagesc(movieSignal(:,:,1))
colormap('gray'); %how the sitmulus was shown
axis equal; axis tight;
xticks([]); yticks([]);

%plot the prediction trace
subplot(2,1,2)
%only for visual comparison we normalize the predicted trace to have the
%same mean and std as the response to the movie.
x=predMemP_y;
m2=mean(memP_AVGMovie);
s2=std(memP_AVGMovie);
s1=std(x);
m1=mean(x);
corrPredMemPred=m2+(x-m1)*(s2/s1);

plot(memP_AVGMovie(STA_output.STA_fnbr:end),'k') %original membrane potential
hold on;
plot(corrPredMemPred,'r')%predicted membrane potential
title(['pCorr:' mat2str(round(pearCorrSqrt,2)) ' Filename: ' STA_output.filename])
xlabel('stimulus frames')
ylabel('voltage (ms)')
legend('original response','predicted response')

%-------------------------------------------------------------------------------
%% 5. add variables into the output
naturalMovie_output.pearCorrSqrt=pearCorrSqrt;
naturalMovie_output.memP_AVGMovie=memP_AVGMovie;
end

function [IndS,IndL]= getIndex(x,valConv,Points)
IndS=find(x<valConv,Points,'last'); %find the last 5 closest smaller values
IndL=find(x>=valConv,Points);
if isempty(IndS)|| length(IndS)<Points %if you are on the left border
    if isempty(IndS) %take value from 1:10
        IndL=find(x>=valConv,Points*2);
        IndS=IndL(1);
        IndL=IndL(2:end);
    else
        IndL=find(x>=valConv,Points*2-length(IndS));
    end
elseif isempty(IndL)|| length(IndL)<Points %if you are on the right border
    if isempty(IndL) %take value from 1:10
        IndS=find(x<valConv,Points*2,'last');
        IndL=IndS(end);
        IndS=IndS(1:end-1);
    else
        IndS=find(x<valConv,Points*2-length(IndL));
    end
end
end

% %make SVD of the filter
% [sp,sigma,temp] = svd(STA_output.STA2D_zoom,'econ');
% sp=sp(:,1);
% imagesc(reshape(sp,size(STA_output.rangeY_zoom,2),size(STA_output.rangeX_zoom,2)));
% temp=temp(:,1);
% %linear filters
% linearFiler1=normSTA(sp);
% linearFiler2=normSTA(temp);
% stim2pred=filter2(linearFiler1,movieSignal2D_zoom,'valid');
% stim2pred=filter2(linearFiler2',stim2pred,'valid');