%this script present the behavior of every participant in homing
%they update positions from path integration
%they combine the headings from Path integration and Landmarks
%they then use updated heading and position to update the home location
%In each cue condition, we simulate N trials
%We also verify this simulation by comparing the predicted homing error SD
%with those predicted by the equations
%at last we compare the predicted SD with the observed SD
%the data set is from the L2/L1=2 condition in Experiment 1

clear all;
close all;

%Observed SD [heading_PI heading_LM homing_PI homing_LM	homing_BOTH	homine_Conflict]	
%[He_SD_PI	He_SD_LM	pi_SD	lm_SD	bo_SD	co_SD
%this data is from L2/L1=2 in Experiment 1
myData=[	
27.798	44.188	12.297	40.153	22.006	18.491	;
26.169	6.477	15.107	27.341	28.252	82.055	;
21.095	23.36	23.3	49.467	27.744	27.266	;
47.03	18.307	28.193	8.358	22.843	22.931	;
25.367	10.047	20.502	11.433	13.247	16.645	;
20.274	8.858	7.194	33.922	29.096	20.147	;
52.366	8.767	20.245	83.356	72.316	61.167	;
31.588	28.837	18.085	37.142	19.8	33.819	;
37.087	29	20.797	10.402	27.196	27.697	;
126.437	13.553	14.779	71.811	62.204	84.657	;
61.946	8.232	39.429	80.778	25.262	27.084	;
39.448	13.652	49.865	70.368	90.474	46.596	;
26.818	2.922	21.105	18.89	11.406	18.627	;
42.703	7.481	25.048	18.722	17.09	19.807	;
10.929	10.622	21.06	13.209	29.08	29.283	;
43.704	23.615	9.882	72.068	121.796	41.61	;
27.187	3.888	5.299	7.764	18.196	42.216	;
30.438	28.349	24.559	29.996	25.934	30.965	;
38.68	24.383	13.316	19.263	34.04	45.851	;
20.058	16.318	15.336	37.32	29.692	37.738	;
18.077	13.942	24.822	42.163	46.141	12.668	;
46.231	20.035	18.274	19.167	24.229	36.15	;
35.947	18.569	5.472	29.817	12.033	8.545	;
11.406	3.808	20.963	11.879	28.967	22.598	;
88.62	32.467	37.287	20.334	16.59	54.841	;
17.476	20.533	14.174	21.117	6.513	37.316	;
25.553	27.183	17.795	12.215	30.082	15.359	;
22.176	34.096	16.436	28.438	2.086	39.475	;
	];


predicted_homing_SD_simulation=zeros(size(myData, 1),4); %save for the predicted SD, will be compared to the observed SD in last four column of myData
predicted_homing_SD_equation=zeros(size(myData, 1),4); %save for the predicted SD, will be compared to the observed SD in last four column of myData


for i = 1: size(myData, 1) 
    % for each subject
%Model
%homing_Error = position_Error - heading_Error see equation 5
%position_Error=(a× heading_Error_PI + ue)  see equation 9
%heading_Error=W_PI × heading_Error_PI+ W_LM× heading_Error_LM see equatino 1
%homing_Error = (a× heading_Error_PI + ue) – (W_PI× heading_Error_PI+ W_LM× heading_Error_LM) see equation 10

%inputs
%He_SD_LM	He_SD_PI; 
%observed SD for heading estimations in Landmark and
%Path-Integration condition from each participant
%heading_Error_PI =normrnd(0, HeadError_PI_SD, N)
%heading_Error_LM =normrnd(0, HeadError_LM_SD, N)

%a and SD of unexplained error
%a      UE Variance, got from the regression line (see equation 9, regression line in Figure 12A and table 8)
%0.79	406.99

%mean optimal weights calcuated from the propostions observed variance of heading errors (PI/(PI+LM)) to predict for the both condition 
%mean observed weight in heading estiamtion in the Conflict condition to predict for the Conflict condition
%optimal_w_LM_HE	obs_w_LM_HE
%0.73415265         0.648131746 weight for heading

%output based on the model
%predicted SD for homing error 

%below is the implementation
% Simulation
% column is in the order of [PI LM Both Conflict]

%observed SD for heading estimations in Landmark and Path-Integration
%conditions
%He_SD_LM	He_SD_PI	
HeadingError_PI_SD= myData(i, 1); %subject i, first column is PI 
HeadingError_LM_SD= myData(i,2);

%a and SD of unexplained error
%a      Varianceue, got from the regression line (see regression line and table)
%0.79	406.99
ue_SD=sqrt(406.99);
a=0.79;

%mean optimal weights (using variance proportion from PI and LM conditions) and mean observed weight for heading errors in the Conflict condition (from group mean)
%optimal_w_LM_HE	obs_w_LM_HE
 optimal_w_LM_HE=0.73415265;       
 obs_w_LM_HE=0.648131746;

% Landmark Weight for heading estimation in conditions of [PI LM Both Conflict]
% column is [PI LM Both Conflict]
W_LM = [0, 1, optimal_w_LM_HE ,  obs_w_LM_HE]; 

W_PI=1-W_LM; % PathIntegration Weight should be [1, 0, optimal_w_PI_HE ,  obs_w_PI_HE]; 

N=100000;

%?PI heading errors in estimating the turning angles which should affect the position errors in all conditions
% column is [PI LM Both Conflict]
% PI heading errors in the second column will time W_PI(2)=0 and not contribute to the combined heading
HeadingError_PI=normrnd(0, HeadingError_PI_SD, [N 4]);

%?LM errors in estimating the heading using the landmark 
% column is [PI LM Both Conflict]
%the first column (for PI) will time W_LM(1)=0 so not contribute to the combined heading
HeadingError_LM=normrnd(0, HeadingError_LM_SD, [N 4]); 
ue=normrnd(0, ue_SD, [N 4]);%unexplained error

% column is [PI LM Both Conflict]
PositionError=a*HeadingError_PI + ue;

% combination of heading errors in all conditions
% column is [PI LM Both Conflict], 
%the first column only contributed from PI heading error as W_LM(1)=0,W_PI(1)=1;  
%the second column only contributed from LM heading error as W_LM(2)=1,W_PI(2)=0; 
% for each row and each column, HeadingError=HeadingError_LM*W_LM +HeadingError_PI * W_PI
HeadingError=HeadingError_LM*diag(W_LM )+HeadingError_PI * diag (W_PI);

%predicted HomingError for each trial and for each condition (N By 4)
HomingError= PositionError -HeadingError;


%%%verfication
% we can calculate homing error SD in all four conditions using equations as well. 
%They should be the same if N in the above simulation is very larger, see
%Table 7B
HomingError_PI_SD= sqrt((a-1)^2 * HeadingError_PI_SD^2 + ue_SD^2);
HomingError_LM_SD= sqrt(a^2*HeadingError_PI_SD^2 +HeadingError_LM_SD^2 + ue_SD^2);
HomingError_Both_SD= sqrt((a-W_PI(3))^2*HeadingError_PI_SD^2 + W_LM(3)^2*HeadingError_LM_SD^2+ ue_SD^2);
HomingError_Conflict_SD= sqrt((a-W_PI(4))^2*HeadingError_PI_SD^2 + W_LM(4)^2*HeadingError_LM_SD^2+ ue_SD^2);

%show the results from simulation and equations
%difference between them should be close to 0

%c=[a;b]
HomingError_sd_simulation=std(HomingError);
HomingError_sd_equation=[HomingError_PI_SD, HomingError_LM_SD,HomingError_Both_SD, HomingError_Conflict_SD];
HomingError_sd_simulation-HomingError_sd_equation

predicted_homing_SD_simulation(i,:)=HomingError_sd_simulation;
predicted_homing_SD_equation(i,:)=HomingError_sd_equation;
end 

%just to plot the predicted from simulation and equations, should fit
%perfectly if N is large
figure 
O=mean(predicted_homing_SD_equation);
P=mean(predicted_homing_SD_simulation);
y = [O', P'];
b = bar(y);
legend(b,'equation','simulation')
set(gca,'XTickLabel',{'PI','LM','Both','Conflict'});



%just to plot the predicted and observed homing SD together, it is similar
%to the result figure of modelling for L2/L1=2 in Experiment 1. Figure 13A
figure
O=mean(myData(:,3:6));
P=mean(predicted_homing_SD_simulation);
y = [O', P'];
b = bar(y);
legend(b,'observed','predicted')
set(gca,'XTickLabel',{'PI','LM','Both','Conflict'});
title('L2/L1=2')


axislim=45;
%xlim([0 axislim])
ylim([0 axislim])