%%  Documentation
%script to analyse MOT experiment data this script loads the experimental
%data into a single structure, fits the weibull cdf to the data,
%computes the best fit model and parameters given the directory the simulations are saved in.
%note, subject 6 was excluded due to having extremely poor task performance
%in one session.

%% set corresponding paths
parentpath = cd(cd('..'));
    datadir=fullfile(parentpath,'data/subdata/matfiles');
    simdir=[fullfile(parentpath,'/data/simulationdata/FE  ')...
        ;fullfile(parentpath,'/data/simulationdata/NEcd')...
        ;fullfile(parentpath,'/data/simulationdata/NEid');...
        fullfile(parentpath,'data/simulationdata/NV  ')];

%% Load subject data and fit weibull to them
cd(datadir) %DATA directory
S=dir('*.mat');
b=[S(:).bytes];
[~,b]=sort(b,'ascend');
b={S(b).name};
E=[];
dist=1.75;

for i =1:length(b)
    tmp=char(b(i));
    if strcmp(tmp(1:2),'S6')==0 %exclude subject 6
    load(char(b(i)));
    fprintf('\n loading file %s \n',tmp)
    sessiondata.sigVcond=flipud(unique(sessiondata.sigV));
    sessiondata.sigXcond=sessiondata.sigX(1);
    sessiondata.sub=sessiondata.subject(1);
    sessiondata.dim= sessiondata.dimension(1);
        for j=1:length(sessiondata.sigVcond)
            
        sessiondata.propCor(j,:)=mean(sessiondata.correct(sessiondata.sigV==sessiondata.sigVcond(j)));
        end
        
        

    
    sessiondata.StimLevels =1./fliplr(m2visAngle(dist,sessiondata.sigVcond'));
    sessiondata.OutOfNum = 30*ones(size(sessiondata.StimLevels));
    sessiondata.NumPos = fliplr(sessiondata.propCor'.*sessiondata.OutOfNum);
    sessiondata.paramsFree = [1 1  0 1]; %estimate threshold, slope, lapse, and fix guess
    
    %search through this grid of parameter values in iterative parameter search
    searchGrid.alpha = 1./[.01:.005:.5]; %threshold
    searchGrid.beta = 10.^[-1:.01:1]; %slope
    searchGrid.gamma = .5; %guess rate fixed
    searchGrid.lambda = [linspace(0,0.05,20)]; %lapse rate
    
    %fit Weibull psychometric curve to data
    [sessiondata.paramsFittedW, LLW] = PAL_PFML_Fit(sessiondata.StimLevels,sessiondata.NumPos,sessiondata.OutOfNum,searchGrid,sessiondata.paramsFree, @PAL_Weibull,'lapseLimits',[0 0.2]);
    sessiondata.weiThresh=sessiondata.paramsFittedW(1);
    sessiondata.weiSlope=sessiondata.paramsFittedW(2);
    sessiondata.weiGuess=sessiondata.paramsFittedW(3);
    sessiondata.weiLapse=sessiondata.paramsFittedW(4);
    x=1./linspace(0.001,0.21,50);
    sessiondata.NumPos=sessiondata.NumPos';
    sessiondata.StimLevels=1./sessiondata.StimLevels';
    
    sessiondata.weiCurve=PAL_Weibull(sessiondata.paramsFittedW,x)'; % creat psychometric curve predictions
    compPoint=0.75; %comparison point of psychometric curves
	% use eq () to compute the subject velocity threshold
    sessiondata.wei80pThresh=(1./(sessiondata.weiThresh *(log((sessiondata.weiLapse+0.5-1)/((sessiondata.weiLapse+compPoint-1))))^(1/sessiondata.weiSlope)));
    E=addstruct(E,sessiondata,'column','force');
    fprintf('psychometric curve fitted for %s\n',tmp)
    end
end
data.subData=E;
fprintf('finished fitting subjects psychometric curves \n ')

%% load in simulation data and compute likelihood
sca=100;
modelInf={'FE','NE-0.90','NE-00','NV'};
modelInf=cellstr(modelInf);
nSubs=sum(length(unique(E.sub)));
cond=unique(E.sigXcond);
nSub=length(unique(E.sub));
E.subCor2D=nan(10,nSub,3);
E.subCor3D=nan(10,nSub,3);
for h=1:length(unique(E.sigXcond))
   data.subData.meanCor2D(:,h)=nanmean(E.propCor(:,E.sigXcond==cond(h)&E.dim==0),2); %calculate average proportion correct
   data.subData.meanCor3D(:,h)=nanmean(E.propCor(:,E.sigXcond==cond(h)&E.dim==1),2);
   data.subData.stdCor2D(:,h)=nanstd(E.propCor(:,E.sigXcond==cond(h)&E.dim==0),0,2);
   data.subData.stdCor3D(:,h)=nanstd(E.propCor(:,E.sigXcond==cond(h)&E.dim==1),0,2);
end

direct=cellstr(simdir);
data.model.parC=nan(60,4);
data.model.parD=nan(60,4);
fprintf('computing model likelihoods \n ')

for r=1:length(direct)
    cd(direct{r})
    
    S=dir('*.mat');
    [C,I]=sort({S.date});
    clear modelAllresp
    clear totalLikeli
    clear LikCor2D
    clear LikCor3D
    for u=1:length(I)
        load(S(I(u)).name);
        data.model.parC(u,r)=C;
        data.model.parD(u,r)=D;
        for j=1:size(data.subData.meanCor2D,1)
            for k=1:size(data.subData.meanCor2D,2)
                if answer2D(k,j)==1
                    answer2D(k,j)=answer2D(k,j)-0.00001;
                end
                if answer3D(k,j)==1
                    answer3D(k,j)=answer3D(k,j)-0.00001;
                end
                LikCor2D(j,k)=log(binopdf(int16(data.subData.meanCor2D(j,k)*(E.OutOfNum(1)*nSub)),(E.OutOfNum(1)*nSub),answer2D(k,j)));
                LikCor3D(j,k)=log(binopdf(int16(data.subData.meanCor3D(j,k)*(E.OutOfNum(1)*nSub)),(E.OutOfNum(1)*nSub),answer3D(k,j)));
            end
        end
        totalLikeli(u)=sum(sum(LikCor2D+LikCor3D));
        modelAllresp(:,:,u)=[answer2D answer3D];
    end
    storeModdelResp{r}=modelAllresp;
    data.model.modelLikeli(r)={totalLikeli};
    [ d, f]=sort(-data.model.modelLikeli{r},'ascend');
    data.model.modelBestResp(:,:,r)=modelAllresp(:,:,f(1));
end
fprintf('finished computing model likelihoods \n ')
%%

%% fit psychometric curve to model data
fprintf('fitting psychometric curve to model data \n ')
for r =1:length(direct)
    for k=1:size(data.subData.meanCor2D,2)
        StimLevels=1./m2visAngle(dist,flipud(unique(E.sigVcond(:,E.sigXcond==cond(k)))));
        OutOfNum = 1000*ones(size(StimLevels)); % we used 1000 simulations so treat this as number of trials
        NumPos2d = data.model.modelBestResp(k,1:10,r);
        NumPos3d = data.model.modelBestResp(k,11:end,r);
        paramsFree = [1 1  0 1]; %estimate threshold, lapse and slope, fix guess

        %fit Weibull
        [paramsFittedW, ~] = PAL_PFML_Fit(StimLevels,NumPos2d'.*OutOfNum,OutOfNum,searchGrid,paramsFree, @PAL_Weibull,'lapseLimits',[0 0.3]);
        data.model.weiThreshMod2D=paramsFittedW(1);
        data.model.weiSlopeMod2D=paramsFittedW(2);
        data.model.weiGuessMod2D=paramsFittedW(3);
        data.model.weiLapseMod2D=paramsFittedW(4);
        x=1./linspace(0.001,0.21,30);
        data.model.weiCurve2D(:,k)=PAL_Weibull(paramsFittedW,x)';
      
        [paramsFittedW, LLW] = PAL_PFML_Fit(StimLevels,NumPos3d'.*OutOfNum,OutOfNum,searchGrid,paramsFree, @PAL_Weibull,'lapseLimits',[0 0.3]);
        data.model.weiThreshMod3D=paramsFittedW(1);
        data.model.weiSlopeMod3D=paramsFittedW(2);
        data.model.weiGuessMod3D=paramsFittedW(3);
        data.model.weiLapseMod3D=paramsFittedW(4);
        data.model.wei80pThreshMod3D(k,r)=real(1./(data.model.weiThreshMod3D *(log((data.model.weiLapseMod3D+0.5-1)/((data.model.weiLapseMod3D+compPoint-1))))^(1/data.model.weiSlopeMod3D)));
        data.model.wei80pThreshMod2D(k,r)=real(1./(data.model.weiThreshMod2D *(log((data.model.weiLapseMod2D+0.5-1)/((data.model.weiLapseMod2D+compPoint-1))))^(1/data.model.weiSlopeMod2D)));
        
        x=1./x;
        
    end
    
    
end
fprintf('finished fitting psychometric curve to model data \n ')
cd(fullfile(parentpath,'/data/subdata/Summary'))
save('motData.mat', 'data')
cd(fullfile(parentpath,'Analysis'))
clear all


