% Main files, describing the essential steps in de model of McGuire and
% Sabes, Nature Neurosciesnce, 2009. 
% The model simulates several data sets, among which the data of Beurze et
% al. J Neurophys, 2006.


% set variance and mu of native signals and prior
SIGMAs = [10 10 10 10 10];    % Sigma's gaze, tv, fv, tp, fp,
MUs    = [-30 30 50 0 20]; % MU's gaze, tv, fv, tp, fp,
SIG_gp = 10; % prior on gaze relative to target

select=0; % single trial (0) or multiple trial (1) version

if select==1, %make MUs matrix for Beurze data
    M=[];
    gaze=[-30:10:30]; tp=[-30:10:30]; fp=[-30:10:30];
    MUs=combvec(gaze,tp,fp)'; MUs(:,4:5)=MUs(:,2:3);
    MUs(:,2)=MUs(:,4)-MUs(:,1);MUs(:,3)=MUs(:,5)-MUs(:,1);
    MUs(find(abs(MUs(:,2))>40 | abs(MUs(:,4)-MUs(:,5))>40 | MUs(:,4)==MUs(:,5)),:)=[];
end
 
x=[-100:.1:100]; % define location space
         
for i=1:size(MUs,1);
                
% native signals

%% representation of gaze 
g=MUs(i,1);  
sig_g=SIGMAs(1); %sig_g=0.1*sig_g+2;
P_g = normpdf(x,g,sig_g);

%% visual representation of target location
tv=MUs(i,2); 
sig_tv=SIGMAs(2);sig_tv=0.1*abs(tv)+3;
P_tv = normpdf(x,tv,sig_tv);

%% visual representation of fingertip location
fv=MUs(i,3);
sig_fv=SIGMAs(3);% sig_fv=0.1*abs(fv)+3;
P_fv = normpdf(x,fv,sig_fv);

%% proprioceptive representation of target location
tp=MUs(i,4); %prop target re body
sig_tp=SIGMAs(4);
P_tp = normpdf(x,tp,sig_tp);

%% proprioceptive representation of fingertip location
fp=MUs(i,5); 
sig_fp=SIGMAs(5);
P_fp = normpdf(x,fp,sig_fp);

%% prior on gaze to tp
P_gp=normpdf(x, tp, SIG_gp); %prior on gaze
P_g_gp=P_g.*P_gp; P_g_gp=P_g_gp/sum(P_g_gp); %gaze biased to target

%% transformed target signals, biased by gaze prior
P_tp_RET=conv(P_tp, P_g_gp(end:-1:1),'same'); %transformed tp
P_tv_BODY=conv(P_tv, P_g_gp,'same'); %transformed tv

%% transformed finger location signals, unbiased 
P_fv_BODY=conv(P_fv, P_g,'same'); 
P_fp_RET=conv(P_fp, P_g(end:-1:1),'same'); 

%%integrate native and transformed signals:
P_finger_RET=P_fv.*P_fp_RET; P_finger_RET=P_finger_RET/sum(P_finger_RET);
P_target_RET=P_tv.*P_tp_RET; P_target_RET=P_target_RET/sum(P_target_RET);

P_finger_BODY=P_fp.*P_fv_BODY;P_finger_BODY=P_finger_BODY/sum(P_finger_BODY);
P_target_BODY=P_tp.*P_tv_BODY;P_target_BODY=P_target_BODY/sum(P_target_BODY);

%% compute displacement vectors
P_MV_RET=conv(P_target_RET, P_finger_RET(end:-1:1),'same'); 
P_MV_BODY=conv(P_target_BODY, P_finger_BODY(end:-1:1),'same'); 
P_MV_INTEG=P_MV_RET.*P_MV_BODY; P_MV_INTEG=P_MV_INTEG/sum(P_MV_INTEG);

%% compute MAP's
MV_RET=x(find(P_MV_RET==max(P_MV_RET)));
MV_BODY=x(find(P_MV_BODY==max(P_MV_BODY)));
MV_INTEG=x(find(P_MV_INTEG==max(P_MV_INTEG)));


M=[M; MUs(i,:) MV_RET MV_BODY MV_INTEG];

end 

close all


%% Plotting routines for multiple trial version
if select==1,
    figure
    subplot(1,2,1); hold on;
    P=[];i=1
    for T=[-30 -20 -10 0 10 20 30];
        MT=M(find(M(:,4)==T),:);
        j=1;
        for TE=[-40:10:40];
            ind=find(MT(:,2)==TE);
            if ~isempty(ind), P(i,j)=mean(MT(ind,4)-MT(ind,5)-MT(ind,7)); else P(i,j)=NaN; end
            j=j+1;
        end
        i=i+1;
    end
    plot(-40:10:40,P,'o-');
    xlabel('Eye Re Target');
    ylabel('Error');
    
    subplot(1,2,2); hold on;
    P=[];i=1
    for T=[-30:10:30];
        MT=M(find(M(:,4)==T),:);
        j=1;
        for TE=[-40:10:40];
            ind=find((MT(:,4)-MT(:,5))==TE);
            if ~isempty(ind), P(i,j)=mean(MT(ind,4)-MT(ind,5)-MT(ind,7)); else P(i,j)=NaN; end
            j=j+1;
        end
        i=i+1;
    end
    
    plot(-40:10:40,P,'o-');
    xlabel('Hand Re Target');
    ylabel('Error');
end