classdef collectiveProj3620220121 < handle
   %version 20220120 added the noise in the simulation that produces the
   %resp location.   %noise to the simulated locations
   %getSiOresp
   %simulationPara has 9 parameters, the last one is the noise %rmse
   %version 20220121, we have one noise sheet inclding all trials(896), which is
   %associated with ID in the split sheets. Across all 5 cv,the same noises are used for
   %each ID,similiar to the observed Oresp that is fixed for each ID.
    properties (SetAccess = private)
        objectTotal=0;
    end
    
    properties (SetAccess = private)
        whichObjects=[0];
    end
    
    properties (SetAccess = private)
        legRatio=0;
    end
    
    
    properties (SetAccess = private)
        VectorCal=collectiveVectorCal;
    end
    properties (SetAccess = private)
        IterationNumber=10;
    end
    properties (SetAccess = private)
        simulationPara=[];
    end
    properties (SetAccess = private)
        %the encoPara from encoding only in Step 3 can be used to calcuate Ppred and hpred
        %that are used to fit the response only in Step 3
        encoParaTable=table;
    end
    properties (SetAccess = private)
        %we need to save noise for each sheet based on the sheet name, so
        %different models will use the same noise to do fitting (3 models so 3 times), and then
        %do verification using the same or different parameters to do
        %validation (3X2)
        noise=[]%containers.Map; %
    end
    methods
        function self=collectiveProj3620220121(ObjTl,ObjID,IterNumber,para)       
            if nargin == 0
                ObjTl = 0;
                ObjID=[0];
                IterNumber=10;
                para=[];
            elseif nargin == 1
                ObjID=[0];
                IterNumber=10;
                para=[];
            elseif nargin == 2 
                IterNumber=10;
                para=[];
            elseif nargin == 3 
                para=[];
            elseif nargin > 4
                print("too much parameters");
            end
            self.objectTotal=ObjTl;
            self.whichObjects=ObjID;
            self.IterationNumber=IterNumber;
            self.simulationPara=para;
            %newpath='/Users/yafeiqi/Downloads/Modelling/Summary/version220121/CircStat2012a/'
            newpath='\\weiminmlab\labmeeting\labmeeting\circular statistics\circular statistic for matlab\CircStat2012a'
            userpath(newpath);
            
            sheetName={};
            encoPara=[];
            t1=table(sheetName,encoPara);
            self.encoParaTable=[self.encoParaTable, t1];%(sheetName, encoPara);
            
        end
        
        function [P,T,O,Oresp,alpha,OT,TP,v_OT,Pprime,hprime,h,Ppred,hpred]=getPathInfo(self,fileName,sheetName,flag,objectTotal,WhichObjects,legRatio)
            
            fTable= readtable(fileName,'Sheet',sheetName);
           
            P=[fTable.testing_positionX,fTable.testing_positionZ];
            T=[fTable.Turning_X,fTable.Turning_Z];              %%%%% Turning point
            Pprime=[fTable.PprimeBIX, fTable.PprimeBIY]; 
            hprime=[fTable.hprimeBI]; 
            
              
            O=[fTable.Origin_X,fTable.Origin_Z];      %%%%% Origin

            Oresp=[fTable.response_O_position_X,fTable.response_O_position_Z]; %% actual response of origin

            for i=1:objectTotal
                columnX=find(ismember(fieldnames(fTable),['obj',num2str(i),'_pos_X']));
                columnZ=find(ismember(fieldnames(fTable),['obj',num2str(i),'_pos_Z']));


                X=fTable{:,columnX};
                Z=fTable{:,columnZ};
                O=[O,X,Z];

                %%response_obj1_X
                columnX=find(ismember(fieldnames(fTable),['response_obj',num2str(i),'_X']));
                columnZ=find(ismember(fieldnames(fTable),['response_obj',num2str(i),'_Z']));


                X=fTable{:,columnX};
                Z=fTable{:,columnZ};
                Oresp=[Oresp,X,Z];
            end
            
            [O, Oresp]=self.selectOOresp(O,Oresp,WhichObjects);
            
            

            alpha=fTable.turnAngle; 
         
            OT=[fTable.OT];     %%%% length of the first leg
            TP=[fTable.TP];   %%%% length of the second leg

            v_OT=self.calculate_vector(O(:,1:2),T); 
            
            v_TP=P-T;
            h=self.VectorCal.getangle_vAtovB([0,1],v_TP);
            
            if strcmp(flag,'resp')
                [Ppred,hpred]=self.getPpredhpred(sheetName,O, OT, TP, v_OT, alpha);% different from Pprime, hprime
            else
                Ppred=Pprime;
                hpred=hprime;
            end
%                 Ppred=Pprime;
%                 hpred=hprime;
            if legRatio
                LegRatio=TP./OT;
                %index=find(LegRatio<=1);
                index=find(LegRatio<=1);
                P=P(index,:);T=T(index,:);O=O(index,:);Oresp=Oresp(index,:);alpha=alpha(index,:);OT=OT(index,:);TP=TP(index,:);v_OT=v_OT(index,:);Pprime=Pprime(index,:);
                hprime=hprime(index,:);h=h(index,:);Ppred=Ppred(index,:);hpred=hpred(index,:);
                
            end 
            
        end
        
        %%%add experiment information to differentiate the same paths in different experiments
        
         function [expno,P,T,O,Oresp,alpha,OT,TP,v_OT,Pprime,hprime,h,Ppred,hpred]=getPathInfoWithExpNo(self,fileName,sheetName,flag,objectTotal,WhichObjects,legRatio)
            
            fTable= readtable(fileName,'Sheet',sheetName);
            expno=fTable.Identity;
            P=[fTable.testing_positionX,fTable.testing_positionZ];
            T=[fTable.Turning_X,fTable.Turning_Z];              %%%%% Turning point
            Pprime=[fTable.PprimeBIX, fTable.PprimeBIY]; 
            hprime=[fTable.hprimeBI]; 
            
              
            O=[fTable.Origin_X,fTable.Origin_Z];      %%%%% Origin

            Oresp=[fTable.response_O_position_X,fTable.response_O_position_Z]; %% actual response of origin

            for i=1:objectTotal
                columnX=find(ismember(fieldnames(fTable),['obj',num2str(i),'_pos_X']));
                columnZ=find(ismember(fieldnames(fTable),['obj',num2str(i),'_pos_Z']));


                X=fTable{:,columnX};
                Z=fTable{:,columnZ};
                O=[O,X,Z];

                %%response_obj1_X
                columnX=find(ismember(fieldnames(fTable),['response_obj',num2str(i),'_X']));
                columnZ=find(ismember(fieldnames(fTable),['response_obj',num2str(i),'_Z']));


                X=fTable{:,columnX};
                Z=fTable{:,columnZ};
                Oresp=[Oresp,X,Z];
            end
            
            [O, Oresp]=self.selectOOresp(O,Oresp,WhichObjects);
            
            

            alpha=fTable.turnAngle; 
         
            OT=[fTable.OT];     %%%% length of the first leg
            TP=[fTable.TP];   %%%% length of the second leg

            v_OT=self.calculate_vector(O(:,1:2),T); 
            
            v_TP=P-T;
            h=self.VectorCal.getangle_vAtovB([0,1],v_TP);
            
            if strcmp(flag,'resp')
                [Ppred,hpred]=self.getPpredhpred(sheetName,O, OT, TP, v_OT, alpha);% different from Pprime, hprime
            else
                Ppred=Pprime;
                hpred=hprime;
            end
%                 Ppred=Pprime;
%                 hpred=hprime;
            if legRatio
                LegRatio=TP./OT;
                %index=find(LegRatio<=1);
                index=find(LegRatio<=1);
                P=P(index,:);T=T(index,:);O=O(index,:);Oresp=Oresp(index,:);alpha=alpha(index,:);OT=OT(index,:);TP=TP(index,:);v_OT=v_OT(index,:);Pprime=Pprime(index,:);
                hprime=hprime(index,:);h=h(index,:);Ppred=Ppred(index,:);hpred=hpred(index,:);
                
            end 
            
        end

        

        
        function  [vector] = calculate_vector(self,A,B)
          vector = B - A;
        end
        
        
        function [rsq,  SSresid, Llog,AIC]=regressionWithCoefficents(self,x,y,coefficients,k)%k is 2 if coefficients are from regression or 0 if coefficients are prior independent of data
        %%%https://www.mathworks.com/help/matlab/data_analysis/linear-regression.html
            p=coefficients;
            yfit=polyval(p,x);
            yresid = y - yfit;
            SSresid = sum(yresid.^2);
            SStotal = (length(y)-1) * var(y);
            rsq = 1 - SSresid/SStotal;
            Llog=self.MLEunivar(yresid);
            AIC=2*(k+1)-2*Llog;
        end
        function Llog=MLEunivar(self,yresid)
            mu=0;
            SSresid = sum(yresid.^2);
            Sigma=SSresid/length(yresid);
            Prob = mvnpdf(yresid,mu,Sigma);  %the cumulative distribution function of the multivariate normal distribution
             %with mean mu and covariance Sigma, evaluated at each row.
            Plog= log(Prob);
            Llog=sum(Plog); %product of all p
        end

        function mdl=regressionLine(self,Angle_TPO1, Angle_TPOresp1)
            figure
            %plot(Angle_TPO1, Angle_TPOresp1,'.');
            scatter(Angle_TPO1, Angle_TPOresp1);
            hold on
            mdl_angle=fitlm(Angle_TPO1, Angle_TPOresp1,'linear');

            x=Angle_TPO1;
            plot(x,x,'r','LineWidth',0.7)
            hold on
            Slope = mdl_angle.Coefficients.Estimate(2);
            Intercept = mdl_angle.Coefficients.Estimate(1);
            r = sqrt(mdl_angle.Rsquared.Ordinary);
            ye=mdl_angle.Coefficients.Estimate(2).*x+mdl_angle.Coefficients.Estimate(1);
            %plot(x,mdl_angle.Coefficients.Estimate(2).*x+mdl_angle.Coefficients.Estimate(1),'Marker','+','LineWidth',1.1)
            plot(x,ye,'Marker','+','LineWidth',1.2)
            %FormStr = sprintf('y = Est()', Est,Rsq.Ordinary);
            if Intercept< 0
                formatSpec = 'y = %.3gx%.2g,r=%.2g'; %y=?x+?,r=?. %.3g denotes a value with three digits
            else
                formatSpec = 'y = %.3gx+%.2g,r=%.2g'; %y=?x+?,r=?. %.3g denotes a value with three digits
            end
            FormStr = sprintf(formatSpec,Slope,Intercept,r);
            
            
            legend({'','Diagonal line',FormStr, 'Location','southeast','FontSize','14'})
            %annotation('textbox', [0.2 0.01 0.4 0.4], 'String', FormStr, 'FitBoxToText','on')
            
            axis equal
            %for inbound path length
%             xlim([0,1.8]);
%             ylim([-1,3.3]);
%             xlabel('Predicted length errors') 
%             ylabel('Mean response length erros') 
%             %for inbound turn angle
            xlim([-10,30]);           %xlim([-10,30]); xlim([-3,59]);
            ylim([-20,60]);           %ylim([-30,60]); ylim([-30,80])
            xlabel('Predicted angle errors') 
            ylabel('Mean response angle erros') 

             mdl=mdl_angle;
        end
        
        
        
        function output=circSD(self,A)
            Aradian = circ_ang2rad(A);
           [ang1,cir1] = circ_std(Aradian);
            output = (circ_rad2ang(cir1));
        end
        function output=circMean(self,A)

            output= circ_rad2ang(circ_mean(circ_ang2rad(A)));

        end


        function [AxUni, AyCmean]=circularMeansForUnqueX(self,Ax,Ay)
            AxUni=unique(Ax);
            AyCmean=[];
            for i=1:numel(AxUni)
                Ay1=Ay(find(Ax==AxUni(i)));
                AyCmean1=self.circMean(Ay1);%mean(Ay1);
                AyCmean=[AyCmean; AyCmean1];
            end
            %AyCmean
        end
        
         function [AxUni, Aymean]=MeansForUnqueX(self,Ax,Ay)
            AxUni=unique(Ax);
            Aymean=[];
            for i=1:numel(AxUni)
                Ay1=Ay(find(Ax==AxUni(i)));
                Aymean1=mean(Ay1);%mean(Ay1);
                Aymean=[Aymean; Aymean1];
            end
            %AyCmean
         end
        
         function [AxUni, Aymean]=MedianForUnqueX(self,Ax,Ay)
            AxUni=unique(Ax);
            Aymean=[];
            for i=1:numel(AxUni)
                Ay1=Ay(find(Ax==AxUni(i)));
                Aymean1=median(Ay1);%mean(Ay1);
                Aymean=[Aymean; Aymean1];
            end
            %AyCmean
        end


        function visualizeAnglePairs(self,Angle1,Angle2)
           
            Angle1=self.VectorCal.RightToLeftAngle(Angle1); Angle2=self.VectorCal.RightToLeftAngle(Angle2);
            fg=figure
            r1=1;
            r2=2;
            self.drawCircle(fg,[0,0], r1);
            self.drawCircle(fg,[0,0], r2);
            x1=r1*cosd(Angle1);
            y1=r1*sind(Angle1);
            x2=r2*cosd(Angle2);
            y2=r2*sind(Angle2);

        
            AxUni=unique(Angle1);
            ColorUni=[1:round(256/numel(AxUni)):(numel(AxUni)-1)*round(256/numel(AxUni))+1];
            ColorUni=ColorUni(randperm(length(ColorUni)));

            for i=1:numel(Angle1)
                col=zeros(1,3);
                p=colormap;
                a=find(AxUni==Angle1(i));
       
                col=p(ColorUni(a),:);
                plot([x1(i),x2(i)],[y1(i),y2(i)],'Color',col);
                hold on
            end

        end
        function drawCircle(self,fg,center, r)
            theta = linspace(-pi,pi);
            xc = r*cos(theta)+center(1);
            yc = r*-sin(theta)+center(2);
            plot(xc,yc);
            hold on
            axis equal
        end 
        
        
        function output=cross_ValidateDisAng(self,fileName,sheet1,sheet2,flag) %flag=h,hprime, or hpred
            trainingSheet=sheet1;
            testingSheet=sheet2;
            output12=self.ValidateDisAng(fileName,trainingSheet,testingSheet,flag);
            trainingSheet=sheet2;
            testingSheet=sheet1;
            output21=self.ValidateDisAng(fileName,trainingSheet,testingSheet,flag);
            output.p=[output12.p;output21.p];
            output.meansqerrorTraining=[output12.train.meansqerror;output21.train.meansqerror];
            output.meansqerrorTesting=[output12.test.meansqerror;output21.test.meansqerror];
            output.partialRsquareTraining=[output12.train.partialRsquare;output21.train.partialRsquare];
            output.partialRsquareTesting=[output12.test.partialRsquare;output21.test.partialRsquare];
            output.LogLTraining=[output12.train.LogL;output21.train.LogL];
            output.LogLTesting=[output12.test.LogL;output21.test.LogL];
            output.LrTraining=[output12.train.Lr;output21.train.Lr];
            output.LrTesting=[output12.test.Lr;output21.test.Lr];
            output.AICTraining=[output12.train.AIC;output21.train.AIC];
            output.AICTesting=[output12.test.AIC;output21.test.AIC];
        end
        function output=ValidateDisAng(self,fileName,trainingSheet,testingSheet,flag)%flag=h,hprime, or hpred

                [p_dis, p_ang]=self.ModelTrainingDisAng(fileName,trainingSheet,flag);
                outputTrain=self.modelTestingDisAng(fileName,trainingSheet,p_dis,p_ang,2,flag);%2 means two parameters estimated in the training data
                outputTest=self.modelTestingDisAng(fileName,testingSheet,p_dis,p_ang,0,flag);

                output.p=[p_dis,p_ang];

                output.train.meansqerror=[outputTrain.baseErrorDis,outputTrain.modelErrorDis,outputTrain.baseErrorAng,outputTrain.modelErrorAng];
                output.train.partialRsquare=[outputTrain.partialRsquareDis,outputTrain.partialRsquareAng];
                output.train.LogL=[outputTrain.baseLogLDis, outputTrain.modelLogLDis,outputTrain.baseLogLAng, outputTrain.modelLogLAng];
                output.train.Lr=[outputTrain.LrDis,outputTrain.LrAng];
                output.train.AIC=[outputTrain.baseAICDis, outputTrain.modelAICDis,outputTrain.baseAICAng, outputTrain.modelAICAng];

                output.test.meansqerror=[outputTest.baseErrorDis,outputTest.modelErrorDis,outputTest.baseErrorAng,outputTest.modelErrorAng];
                output.test.partialRsquare=[outputTest.partialRsquareDis,outputTest.partialRsquareAng];
                output.test.LogL=[outputTest.baseLogLDis, outputTest.modelLogLDis,outputTest.baseLogLAng, outputTest.modelLogLAng];
                output.test.Lr=[outputTest.LrDis,outputTest.LrAng];
                output.test.AIC=[outputTest.baseAICDis, outputTest.modelAICDis,outputTest.baseAICAng, outputTest.modelAICAng];

        end 
       
        function    [O,Oresp,Pg,hg,P,alpha,OT,TP,v_OT,h]=getPathInfoDisAng(self,fileName,sheetName,flag,objectTotal,whichObjects,legRatio)  %flag=h,hprime, or hpred 
            
                fTable= readtable(fileName,'Sheet',sheetName);
     
                P=[fTable.testing_positionX,fTable.testing_positionZ];
                T=[fTable.Turning_X,fTable.Turning_Z];              %%%%% Turning point
               
                h=self.VectorCal.getangle_vAtovB([0,1],P-T);
                O=[fTable.Origin_X,fTable.Origin_Z];      %%%%% Origin
                Oresp=[fTable.response_O_position_X,fTable.response_O_position_Z]; %% actual response of origin
                for i=1:objectTotal
                    columnX=find(ismember(fieldnames(fTable),['obj',num2str(i),'_pos_X']));
                    columnZ=find(ismember(fieldnames(fTable),['obj',num2str(i),'_pos_Z']));


                    X=fTable{:,columnX};
                    Z=fTable{:,columnZ};
                    O=[O,X,Z];

                    %%response_obj1_X
                    columnX=find(ismember(fieldnames(fTable),['response_obj',num2str(i),'_X']));
                    columnZ=find(ismember(fieldnames(fTable),['response_obj',num2str(i),'_Z']));


                    X=fTable{:,columnX};
                    Z=fTable{:,columnZ};
                    Oresp=[Oresp,X,Z];
                end
                
                [O, Oresp]=selectOOresp(self,O,Oresp,whichObjects);

                alpha=fTable.turnAngle; 
                %PathLength=[0.9,1.8];
                OT=[fTable.OT];     %%%% length of the first leg
                TP=[fTable.TP];   %%%% length of the second leg
                
                %proj36=collectiveProj36;
                v_OT= self.calculate_vector(O(:,1:2),T); 
                if strcmp(flag,'hprime')   
                   Pprime=[fTable.PprimeBIX, fTable.PprimeBIY]; 
                   hprime=[fTable.hprimeBI]; 
                   Pg=Pprime;
                   hg=hprime;
                elseif strcmp(flag,'hpred')
                   
                   [Ppred,hpred]=self.getPpredhpred(sheetName,O, OT, TP, v_OT, alpha);
             
                   Pg=Ppred;
                   hg=hpred;
                else
                   Pg=P;
                   hg=h;
                   columnh=find(ismember(fieldnames(fTable),'h'));
                    if ~columnh
                        fTable=[fTable table(h)];
                        writetable(fTable,fileName,'Sheet',sheetName);
                    end
                    %fTable1= [fTable table(hpred,PpredX,PpredY)];
                    
                end
                
                if legRatio
                    LegRatio=TP./OT;
                    %index=find(LegRatio<=1);
                    index=find(LegRatio<=1);         
                    
                    O=O(index,:);
                    Oresp=Oresp(index,:);
                    Pg=Pg(index,:);
                    hg=hg(index,:);
                    P=P(index,:);
                    alpha=alpha(index,:);
                    OT=OT(index,:);
                    TP=TP(index,:);
                    v_OT=v_OT(index,:);
                    h=h(index,:);
                end
                
        end 
                
       

        function [p_dis, p_ang,mdl_dis,mdl_ang,dx,dy,Ax,Ay]=ModelTrainingDisAng(self,fileName,sheetName,flag)
            [O,Oresp,Ppred,hpred,P,alpha,OT,TP,v_OT,h]=self.getPathInfoDisAng(fileName,sheetName,flag,self.objectTotal,self.whichObjects,self.legRatio);
            
              if ~isempty(self.simulationPara) %do simulation instead  
                  
                 Oresp= self.getSiOresp(P,O,alpha,OT,TP,v_OT,self.simulationPara,sheetName,self.whichObjects,fileName);
                
              end
            
            [p_dis, p_ang, mdl_dis,mdl_ang,dx,dy,Ax,Ay]=self.modelfitDisAng(O,P,h,Oresp,Ppred,hpred); 
        end
        function [p_dis, p_ang, mdl_dis,mdl_ang,dx,dy,Ax,Ay]=modelfitDisAng(self,O,P,h,Oresp,Pprime,hprime)

                   [dx,dy,Ax,Ay]=dis_angPredResphomingV(self,O,P,h,Oresp,Pprime,hprime);
                    mdl_dis=self.regressionLine(dx, dy);
                    xlabel('Correct Homebound Length (m)','FontSize',12)
                    ylabel('Observed Homebound Length (m)','FontSize',12)
                    legend({'','Diagonal line', 'y = 0.633x + 2.201, r = .243'})
                    xlim([-3,9])
                    ylim([0,20])
                    p_dis=polyfit(dx, dy,1);
                    mdl_ang=self.regressionLine(Ax, Ay);
                    xlabel('Correct Turn Angle (º)','FontSize',12)
                    ylabel('Observed Turn Angle (º)','FontSize',12)
                    legend({'','Diagonal line', 'y = 0.864x + 28.257, r = .632'})
                    xlim([-50 310])
                    p_ang=polyfit(Ax, Ay,1);
                  
        end
        
        function output=modelTestingDisAng(self,fileName,sheetName,p_dis,p_ang,k,flag)
               
                [O,Oresp,Ppred,hpred,P,alpha,OT,TP,v_OT,h]=self.getPathInfoDisAng(fileName,sheetName,flag,self.objectTotal,self.whichObjects,self.legRatio);
                 if ~isempty(self.simulationPara) %do simulation instead  
                      Oresp= self.getSiOresp(P,O,alpha,OT,TP,v_OT,self.simulationPara,sheetName,self.whichObjects,fileName);
                 end
                output=self.modelValidateDisAng(O,P,h,Oresp,Ppred,hpred, p_dis,p_ang,k);
        end
        function output = modelValidateDisAng(self,O,P,h,Oresp,Pprime,hprime,p_dis,p_ang,k)
            
            [dx,dy,Ax,Ay]=dis_angPredResphomingV(self,O,P,h,Oresp,Pprime,hprime);
            
            [rsqModel_dis, SSresidModel_dis,LogLModel_dis,AICmodel_dis]=self.regressionWithCoefficents(dx, dy,p_dis,k);
            [rsqBase_dis,  SSresidBase_dis,LogLBase_dis,AICBase_dis]=self.regressionWithCoefficents(dx, dy,[1 0],0);

            
            [rsqModel_ang, SSresidModel_ang,LogLModel_ang,AICmodel_ang]=self.regressionWithCoefficents(Ax, Ay,p_ang,k);
            [rsqBase_ang,  SSresidBase_ang,LogLBase_ang,AICBase_ang]=self.regressionWithCoefficents(Ax, Ay,[1 0],0);
            n=length(dx);
          
            
            output.modelErrorDis =SSresidModel_dis/n;
            output.baseErrorDis=SSresidBase_dis/n;
            output.modelLogLDis=LogLModel_dis;
            output.baseLogLDis=LogLBase_dis;
            output.modelAICDis=AICmodel_dis;
            output.baseAICDis=AICBase_dis;
            if SSresidModel_ang/n>0.00001
                output.partialRsquareDis=1-SSresidModel_dis/SSresidBase_dis;
                output.LrDis=exp(output.modelLogLDis-output.baseLogLDis);
            else
                output.partialRsquareDis=0;
                output.LrDis=1;
            end
           
            output.modelErrorAng =SSresidModel_ang/n;
            output.baseErrorAng =SSresidBase_ang/n;
            output.modelLogLAng =LogLModel_ang;
            output.baseLogLAng =LogLBase_ang;
            output.modelAICAng =AICmodel_ang;
            output.baseAICAng =AICBase_ang;
            if SSresidModel_ang/n>0.00001
                output.partialRsquareAng =1-SSresidModel_ang/SSresidBase_ang;
                output.LrAng =exp(output.modelLogLAng-output.baseLogLAng);
            else     
                output.partialRsquareAng =0;
                output.LrAng =1;
            end
            
        end
        
        function output=ValidateLocation(self,fileName, trainingSheet,testingSheet,flag)%'all' 'enco' 'resp'
    
            para1=[1,    0,    1,   0];
            para2=[1,    0,    1,   0];
            para3=[1,    0,    1,   0,   1,    0,    1,   0];
            x1=0;x2=0;x3=0;
            if strcmp(flag,'all') 
                [x1]=self.ModelTrainingLocation(fileName,trainingSheet,para1, 1,flag);
                outputTrain1=self.ModelTestingLocation(fileName,trainingSheet,x1,1,flag);
                outputTest1=self.ModelTestingLocation(fileName,testingSheet,x1,1,flag);
                [x2]=self.ModelTrainingLocation(fileName,trainingSheet,para2, 2,flag);
                outputTrain2=self.ModelTestingLocation(fileName,trainingSheet,x2,2,flag);
                outputTest2=self.ModelTestingLocation(fileName,testingSheet,x2,2,flag);
                para3=[x1,x2];
                [x3]=self.ModelTrainingLocation(fileName,trainingSheet,para3, 3,flag);
                outputTrain3=self.ModelTestingLocation(fileName,trainingSheet,x3,3,flag);
                outputTest3=self.ModelTestingLocation(fileName,testingSheet,x3,3,flag);
            elseif strcmp(flag,'enco') 
                [x1]=self.ModelTrainingLocation(fileName,trainingSheet,para1, 1,flag);
                
                
                self.encoParaTable=[self.encoParaTable; {testingSheet,x1}];
                
                
                outputTrain1=self.ModelTestingLocation(fileName,trainingSheet,x1,1,flag);
                outputTest1=self.ModelTestingLocation(fileName,testingSheet,x1,1,flag);

                outputTrain2=outputTrain1;
                outputTest2=outputTest1;
                outputTrain3=outputTrain1;
                outputTest3=outputTest1;
            else
                [x2]=self.ModelTrainingLocation(fileName,trainingSheet,para2, 2,flag);
                outputTrain2=self.ModelTestingLocation(fileName,trainingSheet,x2,2,flag);
                outputTest2=self.ModelTestingLocation(fileName,testingSheet,x2,2,flag);

                outputTrain1=outputTrain2;
                outputTest1=outputTest2;
                outputTrain3=outputTrain2;
                outputTest3=outputTest2;
            end
            output.x1=x1;
            output.x2=x2;
            output.x3=x3;
           
            output.train.meansqerror=[outputTrain1.baseError,outputTrain1.modelError,outputTrain2.modelError,outputTrain3.modelError];
            output.train.partialRsquare=[0, outputTrain1.partialRsquare,outputTrain2.partialRsquare,outputTrain3.partialRsquare];
            output.train.LogL=[outputTrain1.baseLogL, outputTrain1.modelLogL,outputTrain2.modelLogL,outputTrain3.modelLogL ];
            output.train.Lr=[1, outputTrain1.Lr,outputTrain2.Lr,outputTrain3.Lr];

            output.test.meansqerror=[outputTest1.baseError,outputTest1.modelError,outputTest2.modelError,outputTest3.modelError];
            output.test.partialRsquare=[0, outputTest1.partialRsquare,outputTest2.partialRsquare,outputTest3.partialRsquare];
            output.test.LogL=[outputTest1.baseLogL, outputTest1.modelLogL,outputTest2.modelLogL,outputTest3.modelLogL ];
            output.test.Lr=[1, outputTest1.Lr,outputTest2.Lr,outputTest3.Lr];
          
        end 

        function output=cross_ValidateLocation(self,fileName,sheet1,sheet2,flag)%flag='enco', 'resp', 'all'

            trainingSheet=sheet1;
            testingSheet=sheet2;

            output12=self.ValidateLocation(fileName,trainingSheet,testingSheet,flag);
            trainingSheet=sheet2;
            testingSheet=sheet1;
            output21=self.ValidateLocation(fileName,trainingSheet,testingSheet,flag);
            output.x1=[output12.x1;output21.x1];
            output.x2=[output12.x2;output21.x2];
            output.x3=[output12.x3;output21.x3];
            output.meansqerrorTraining=[output12.train.meansqerror;output21.train.meansqerror];
            output.meansqerrorTesting=[output12.test.meansqerror;output21.test.meansqerror];
            output.partialRsquareTraining=[output12.train.partialRsquare;output21.train.partialRsquare];
            output.partialRsquareTesting=[output12.test.partialRsquare;output21.test.partialRsquare];
            output.LogLTraining=[output12.train.LogL;output21.train.LogL];
            output.LogLTesting=[output12.test.LogL;output21.test.LogL];
            output.LrTraining=[output12.train.Lr;output21.train.Lr];
            output.LrTesting=[output12.test.Lr;output21.test.Lr];
            
        end

        function [x]=ModelTrainingLocation(self,fileName,sheetName,para, whichmodel,flag)

            %VectorCal=collectiveVectorCal;
            [P,T,O,Oresp,alpha,OT,TP,v_OT,Pprime,hprime,h,Ppred,hpred]=self.getPathInfo(fileName,sheetName,flag,self.objectTotal,self.whichObjects,self.legRatio);
            
            
            if ~isempty(self.simulationPara) %do simulation instead
               
                Oresp= self.getSiOresp(P,O,alpha,OT,TP,v_OT,self.simulationPara,sheetName,self.whichObjects,fileName);
                
                
                
                [Pprime,hprime]=self.getPprimehprimeBi(fileName,sheetName,self.objectTotal); 
              
               
            end
            %using predicted O as the response O

            %%can I select only the paths with ratio no larger than 1
         
            LegRatio=TP./OT;
            %index=find(LegRatio<=1);
            index=find(LegRatio<=1 | LegRatio>1);
            [x]=self.modelfitLocation(para,O(index,:),P(index,:),OT(index,:),TP(index,:),Oresp(index,:),v_OT(index,:),alpha(index),Pprime(index,:),hprime(index),h(index),Ppred(index,:),hpred(index),whichmodel,flag);
            
            %[x, Base_meansqerror, minfval,fval]=modelfit(para,O,P,OT,TP,Oresp,v_OT,alpha,whichmodel);


        end
        
        
        
        function [x, minfval,fval]=modelfitLocation(self,para,O,P,OT,TP,Oresp,v_OT,alpha,Pprime,hprime,h,Ppred,hpred, whichmodel,flag)
                options = optimset('MaxFunEvals',10e5);

                if whichmodel==3
                    inita1=self.getInitial1(self.IterationNumber,[0 10],[-3 30],[0 10],[-10 90],2);%intia=getInitial(row,limR,limxyz,limS)
                else 
                    inita1=self.getInitial1(self.IterationNumber,[0 10],[-3 30],[0 10],[-10 90],1);%intia=getInitial(row,limR,limxyz,limS)
                end
                inita=para;
                inita=[inita;inita1];
                minfval=10000000;
                %x%fvalall=[]
                xall=[];
                    %ran=getInita() %yaw pitch roll s x y z 
                %minx=1000;
                for i=1:numel(inita(:,1))
                        [x,fval(i)] = fminsearch (@(x)self.myfunLocation(x,O,P,OT,TP,Oresp,v_OT,alpha,Pprime,hprime,h,Ppred,hpred,whichmodel,flag),inita(i,:), options);
                        if fval(i)<minfval
                            minfval=fval(i);
                            minx=x;
                        end
                end
                x=minx;
                   
        end
    
        function f = myfunLocation(self,x,O,P,OT,TP,Oresp,v_OT,alpha,Pprime,hprime,h,Ppred,hpred,whichModel,flag)   
        %       para=[beta1,beta2,beta3,beta4];
             %whichModel=1, encoding only, 2, response only, 3 combined
            
             if whichModel==1
                 a=x(1);
                 b=x(2);
                 c=x(3);
                 d=x(4);
                 p=1;
                 q=0;
                 m=1;
                 n=0;
             elseif whichModel==2
                 a=1;
                 b=0;
                 c=1;
                 d=0;
                 p=x(1);
                 q=x(2);
                 m=x(3);
                 n=x(4);
             else
                 a=x(1);
                 b=x(2);
                 c=x(3);
                 d=x(4);
                 p=x(5);
                 q=x(6);
                 m=x(7);
                 n=x(8);
             end
             minSVal=0;  %minimal possible value for slopes
             maxSVal=10;  %maximal possible value for slopes
             minIValdis=-3;  %minimal possible value for intercept of Distance function

             minIValang=-10;
             maxIValang=90;
             %[0 10],[-3 30],[0 10],[-10 90]

             error=[];
             for i=1:numel(O(1,:))/2
                 index=find(~isnan(O(:,2*i-1)));
                 
                 if strcmp(flag,'all')  
                     Erem_O = self.predict_Location(a,b, c,d, p,q,m,n,O(index,1:2),O(index,2*i-1:2*i),P(index,:), OT(index), TP(index), v_OT(index,:), alpha(index));
                     error =[error;(Oresp(index,2*i-1)-Erem_O(index,1)).^2+(Oresp(index,2*i)-Erem_O(index,2)).^2];
                 elseif strcmp(flag,'enco') 
                    [Tenc,Ppred, EncodingDist_TP,EncodingAng_alpha,henc] = self.predict_P(a,b, c,d, O(index,1:2), OT(index), TP(index), v_OT(index,:), alpha(index));
                    %Ppred=Penc;
                    error =[error;(Pprime(index,2*i-1)-Ppred(index,1)).^2+(Pprime(index,2*i)-Ppred(index,2)).^2];
                 else
                     Erem_O = self.predict_O(p,q,m,n,O(index,2*i-1:2*i),P(index,:), h(index),Ppred(index,:),hpred(index));
                     error =[error;(Oresp(index,2*i-1)-Erem_O(index,1)).^2+(Oresp(index,2*i)-Erem_O(index,2)).^2];           
                 end
             end
        
              penaltyValue=10;
              if a<minSVal || c<minSVal || a> maxSVal || c> maxSVal ||p<minSVal || m<minSVal || p> maxSVal || m> maxSVal
                  penalty=penaltyValue;
              elseif  d<minIValang || d> maxIValang|| n<minIValang || n> maxIValang % for angle intercept smaller than 360
                  penalty=penaltyValue;
              elseif  b<minIValdis || q<minIValdis     % for distance intercept smaller than -1.8(i.e.,0-0.9*2)
                  penalty=penaltyValue;
              else
                  penalty=0;
              end    

              f=mean(error)+penalty;
              RMSE = sqrt(mean(error));


        end
        function [output]=ModelTestingLocation(self,fileName,sheetName,para, whichmodel,flag)
                [P,T,O,Oresp,alpha,OT,TP,v_OT,Pprime,hprime,h,Ppred,hpred]=self.getPathInfo(fileName,sheetName,flag,self.objectTotal,self.whichObjects,self.legRatio);
                if ~isempty(self.simulationPara) %do simulation instead
                    Oresp= self.getSiOresp(P,O,alpha,OT,TP,v_OT,self.simulationPara,sheetName,self.whichObjects,fileName); 
                    [Pprime,hprime]=self.getPprimehprimeBi(fileName,sheetName,self.objectTotal); 
                end
                [output]= self.modelValidateLocation(para,O,P,OT,TP,Oresp,v_OT,alpha,Pprime,hprime,h,Ppred,hpred,whichmodel,flag);
        end
        
        function [OAll,OrespAll,Erem_OAll]=getORPAll(self,x,O,P,OT,TP,Oresp,v_OT,alpha,Pprime,hprime,h,Ppred,hpred,whichModel,flag) % append below vertically
            
              %whichModel=1, encoding only, 2, response only, 3 combined
             if whichModel==1
                 a=x(1);
                 b=x(2);
                 c=x(3);
                 d=x(4);
                 p=1;
                 q=0;
                 m=1;
                 n=0;
             elseif whichModel==2
                 a=1;
                 b=0;
                 c=1;
                 d=0;
                 p=x(1);
                 q=x(2);
                 m=x(3);
                 n=x(4);
             else
                 a=x(1);
                 b=x(2);
                 c=x(3);
                 d=x(4);
                 p=x(5);
                 q=x(6);
                 m=x(7);
                 n=x(8);
             end

              OrespAll=[];
              Erem_OAll=[];
              OAll=[];

              %modelLogL=0;
              %baseLogL=0;
             for i=1:numel(O(1,:))/2
                 index=find(~isnan(O(:,2*i-1)));
                 if strcmp(flag,'all')  %step 2 always uses 'all'
                    
                     Erem_O = self.predict_Location(a,b, c,d, p,q,m,n,O(index,1:2),O(index,2*i-1:2*i),P(index,:), OT(index), TP(index), v_OT(index,:), alpha(index));

                     Erem_OAll=[Erem_OAll;Erem_O];
                     OrespAll=[OrespAll;[Oresp(index,2*i-1),Oresp(index,2*i)]];
                     OAll=[OAll;[O(index,2*i-1),O(index,2*i)]];                 
                 elseif strcmp(flag,'enco') %only for step 3
                    [Tenc,Penc, EncodingDist_TP,EncodingAng_alpha,henc] = self.predict_P(a,b, c,d, O(index,1:2), OT(index), TP(index), v_OT(index,:), alpha(index));
                    Erem_O =Penc;
                    Erem_OAll=[Erem_OAll;Erem_O];
                    OrespAll=[OrespAll;[Pprime(index,2*i-1),Pprime(index,2*i)]];
                    OAll=[OAll;[P(index,1),P(index,2)]]; 
                    %error =[error;(Pprime(index,2*i-1)-Ppred(index,1)).^2+(Pprime(index,2*i)-Ppred(index,2)).^2];
                 else %only for step 3
                     Erem_O = self.predict_O(p,q,m,n,O(index,2*i-1:2*i),P(index,:), h(index),Ppred(index,:),hpred(index));
                     Erem_OAll=[Erem_OAll;Erem_O];
                     OrespAll=[OrespAll;[Oresp(index,2*i-1),Oresp(index,2*i)]];
                     %base o should be affected by pe he already
                     Base_O=self.predict_O(1,0,1,0,O(index,2*i-1:2*i),P(index,:), h(index),Ppred(index,:),hpred(index));                 %OAll=[OAll;[O(index,2*i-1),O(index,2*i)]]; 
                     OAll=[OAll;Base_O];
                 end

             end
        end
        
       
        
        function [output] = modelValidateLocation(self,x,O,P,OT,TP,Oresp,v_OT,alpha,Pprime,hprime,h,Ppred,hpred,whichModel,flag)   
%         %       para=[beta1,beta2,beta3,beta4];
%              %whichModel=1, encoding only, 2, response only, 3 combined
%              if whichModel==1
%                  a=x(1);
%                  b=x(2);
%                  c=x(3);
%                  d=x(4);
%                  p=1;
%                  q=0;
%                  m=1;
%                  n=0;
%              elseif whichModel==2
%                  a=1;
%                  b=0;
%                  c=1;
%                  d=0;
%                  p=x(1);
%                  q=x(2);
%                  m=x(3);
%                  n=x(4);
%              else
%                  a=x(1);
%                  b=x(2);
%                  c=x(3);
%                  d=x(4);
%                  p=x(5);
%                  q=x(6);
%                  m=x(7);
%                  n=x(8);
%              end
% 
%               OrespAll=[];
%               Erem_OAll=[];
%               OAll=[];
% 
%               %modelLogL=0;
%               %baseLogL=0;
%              for i=1:numel(O(1,:))/2
%                  index=find(~isnan(O(:,2*i-1)));
%                  if strcmp(flag,'all')  %step 2 always uses 'all'
%                     
%                      Erem_O = self.predict_Location(a,b, c,d, p,q,m,n,O(index,1:2),O(index,2*i-1:2*i),P(index,:), OT(index), TP(index), v_OT(index,:), alpha(index));
% 
%                      Erem_OAll=[Erem_OAll;Erem_O];
%                      OrespAll=[OrespAll;[Oresp(index,2*i-1),Oresp(index,2*i)]];
%                      OAll=[OAll;[O(index,2*i-1),O(index,2*i)]];                 
%                  elseif strcmp(flag,'enco') %only for step 3
%                     [Tenc,Penc, EncodingDist_TP,EncodingAng_alpha,henc] = self.predict_P(a,b, c,d, O(index,1:2), OT(index), TP(index), v_OT(index,:), alpha(index));
%                     Erem_O =Penc;
%                     Erem_OAll=[Erem_OAll;Erem_O];
%                     OrespAll=[OrespAll;[Pprime(index,2*i-1),Pprime(index,2*i)]];
%                     OAll=[OAll;[P(index,1),P(index,2)]]; 
%                     %error =[error;(Pprime(index,2*i-1)-Ppred(index,1)).^2+(Pprime(index,2*i)-Ppred(index,2)).^2];
%                  else %only for step 3
%                      Erem_O = self.predict_O(p,q,m,n,O(index,2*i-1:2*i),P(index,:), h(index),Ppred(index,:),hpred(index));
%                      Erem_OAll=[Erem_OAll;Erem_O];
%                      OrespAll=[OrespAll;[Oresp(index,2*i-1),Oresp(index,2*i)]];
%                      %base o should be affected by pe he already
%                      Base_O=self.predict_O(1,0,1,0,O(index,2*i-1:2*i),P(index,:), h(index),Ppred(index,:),hpred(index));                 %OAll=[OAll;[O(index,2*i-1),O(index,2*i)]]; 
%                      OAll=[OAll;Base_O];
%                  end
% 
%              end
            [OAll,OrespAll,Erem_OAll]=getORPAll(self,x,O,P,OT,TP,Oresp,v_OT,alpha,Pprime,hprime,h,Ppred,hpred,whichModel,flag); 

             output.modelError =mean( (OrespAll(:,1)-Erem_OAll(:,1)).^2+(OrespAll(:,2)-Erem_OAll(:,2)).^2);
             output.baseError =mean( (OrespAll(:,1)-OAll(:,1)).^2+(OrespAll(:,2)-OAll(:,2)).^2);
             if output.baseError>0.001 %baseError~=0
                 output.partialRsquare=1-output.modelError /output.baseError;
                 output.modelLogL=self.MaxLikelihoodLocation(Erem_OAll,OrespAll);
                 output.baseLogL=self.MaxLikelihoodLocation(OAll,OrespAll);
                 output.Lr=exp(output.modelLogL-output.baseLogL);
             else
                 output.partialRsquare=0;
                 output.modelLogL=self.MaxLikelihoodLocation(Erem_OAll,OrespAll);
                 output.baseLogL=self.MaxLikelihoodLocation(OAll,OrespAll);
                 output.Lr=1;
             end 
         
        end
        function Erem_O = predict_Location(self,a,b, c,d, p,q,m,n,O,Obj,P,OT, TP, v_OT, alpha)%
            
             [Te,Pe, EncodingDist_TP,EncodingAng_alpha,he]=self.encodingErrorLocation(a,b, c,d, O,OT, TP, v_OT, alpha);
           
             Predict_v1=self.responseErrorLocation(p,q,m,n, he, Pe, Obj);
             %VectorCal=collectiveVectorCal;
             Predict_v=self.VectorCal.RotateVecArrayLeft(alpha-EncodingAng_alpha,Predict_v1);
             Erem_O =  Predict_v + P;  %%% predicted pointing O'
        end 
        
        function [Te,Pe,EncodingDist_TP,EncodingAng_alpha, he]=encodingErrorLocation(self,a,b,c,d, O, OT, TP, v_OT, alpha)
          EncodingDist_OT= a*OT+b;      %the encoded distance of the first leg
          EncodingDist_TP = a*TP+b;      %the encoded distance of the second leg 

          v_OTe = v_OT.* repmat(EncodingDist_OT./OT,1,2); % the encoded vector of OT: OT'
          v_OTex = v_OT .* repmat(EncodingDist_TP ./ OT,1,2);   %%%%the extension of the vectorOT'

          %EncodingAng_alpha(i) = c*alpha(i)+d;            %the encoded turning angle
          EncodingAng_alpha= c*abs(alpha)+d; 
    %       if alpha<0
    %         EncodingAng_alpha=-EncodingAng_alpha;
    %       end
          EncodingAng_alpha=EncodingAng_alpha.*(alpha./abs(alpha));
          theta = -EncodingAng_alpha;
          %R1 = [cosd(theta) -sind(theta); sind(theta) cosd(theta)];  %%%% rotation matrix
          %v_TPex = (R1 * (v_OTex'))';     %%%% rotate the extension of the vectorOT' to represent the vectorTP_encoded
          %VectorCal=collectiveVectorCal;
          %v_TPex =VectorCal.RotateVecArray(theta,v_OTex);
          v_TPex =[cosd(theta).*v_OTex(:,1)-sind(theta).*v_OTex(:,2), sind(theta).*v_OTex(:,1)+cosd(theta).*v_OTex(:,2)];
          Te= v_OTe + O(:,1:2);     %%% the encoded turning point
          Pe = v_TPex+Te;      %%%%the coordinates of the P_encoded from the encoding model
            v_TePe=Pe-Te;
            he=self.VectorCal.getangle_vAtovB([0,1],v_TePe);

        end
        
        function Predict_v=responseErrorLocation(self,p,q,m,n,he, Pe,Obj)
          Ideal_Dis = sqrt((Pe(:,1)-Obj(:,1)).^2+(Pe(:,2)-Obj(:,2)).^2);

          ExecuteDist_PO= p * Ideal_Dis +q;  %the executed distance of the homebound path

          v_TePe= [sind(he),cosd(he)];    %the encoded heading h'
          v_PeObj = Obj-Pe;   %%% required pointing vector(PeObj) based on the encoding-error model
          %Ideal_Ang(i) = atan2d(v_TePe(i,1).*v_PeO(i,2)-v_TePe(i,2).*v_PeO(i,1),v_TePe(i,1).*v_PeO(i,1)+v_TePe(i,2).*v_PeO(i,2));
          %VectorCal=collectiveVectorCal;
          Ideal_Ang =self.VectorCal.getangle_vAtovB(v_TePe,v_PeObj);
        
          ExecuteAng_alpha = m*abs(Ideal_Ang)+n; %linear function applied to magnitude            
          ExecuteAng_alpha=ExecuteAng_alpha.*(Ideal_Ang./abs(Ideal_Ang)); %assign the sign
          
          v_Hprime= v_TePe.* repmat(ExecuteDist_PO,1,2);      %the v_H' with the executed length

          beta=ExecuteAng_alpha;
         
          Predict_v=self.VectorCal.RotateVecArrayLeft(beta,v_Hprime);  %
        end
        
        function Llog=MaxLikelihoodLocation(self,Enc_Oprime,Oresp)
            d= Enc_Oprime-Oresp; % d is n by 2 matrix
            d_squared=sum(d.^2,'all')/length(d(:,1));
            mu=[0,0]; %should 0,0
            dx=d(:,1);
            dy=d(:,2);
            %r=corr(dx,dy);
            vdx2=mean(dx.^2);
            vdy2=mean(dy.^2);
            %vdxy=r*sqrt(vdx2*vdy2);
            vdxy= mean(dx.*dy);
            Sigma=[vdx2,vdxy;vdxy,vdy2];   %Covariance matrix of multivariate normal distribution
            %L= (1/det(Sigma)).^(length(alpha)/2);

            Prob = mvnpdf(d,mu,Sigma);  %the cumulative distribution function of the multivariate normal distribution
             %with mean mu and covariance Sigma, evaluated at each row.
            Plog= log(Prob);
            Llog=sum(Plog); %product of all p
            %we can calcuate Maxl using Sigma, the within level C=2 (dx,dy), and
            %independent (between) cases, n=length(dx)
            c=2; 
            n=length(dx);
            Llog1=(-c*n/2)*log(2*pi())+(-c*n/2)-n/2*log(det(Sigma));
        end
        
        function intia=getInitial1(self,row,limS_d,limIntercept_d,limS_a,limIntercept_a,numberEach)%limS_d=[0,5], limIntercept_d=[-15,15],limS_a=[0,5], limIntercept_a=[-90,90]
                numxyz=row;
                %limxyz=5

                a = 50;
                b = 100;
                r = (b-a).*rand(1000,1) + a;

                S_d= (limS_d(2)-limS_d(1)).*rand(row,numberEach)+limS_d(1);
                intercept_d=(limIntercept_d(2)-limIntercept_d(1)).*rand(row,numberEach)+limIntercept_d(1);
                S_a= (limS_a(2)-limS_a(1)).*rand(row,numberEach)+limS_a(1);
                intercept_a=(limIntercept_a(2)-limIntercept_a(1)).*rand(row,numberEach)+limIntercept_a(1);
                %-limS + 2*limS*rand(row,numberEach);
                %limR=pi
                %intercept= -limIntercept + 2*limIntercept*rand(row,numberEach);
                %limS=3
                intia=[];
                for i=1:numberEach
                    intia=[intia, S_d(:,i),intercept_d(:,i),S_a(:,i),intercept_a(:,i)];
                end
        end
        
        function savePprehpre(fileName,sheetName,para)
            [P,T,O,Oresp,alpha,OT,TP,v_OT,Pprime,hprime,h,Ppred,hpred]=self.getPathInfo(fileName,sheetName,flag,self.objectTotal,self.whichObjects,self.legRatio);
            [Tenc,Penc, EncodingDist_TP,EncodingAng_alpha,henc] = self.predict_P(para(1),para(2), para(3),para(4), O, OT, TP, v_OT, alpha);
            PpredX=Penc(:,1);PpredY=Penc(:,2);
            fTable= readtable(fileName,'Sheet',sheetName);
            hpred=henc;
            %vectorCal=collectiveVectorCal;
            hpred=self.VectorCal.ChangeAngleTo_negative180to180(hpred);

            columnhpred=find(ismember(fieldnames(fTable),'hpred'));
            if columnhpred
                fTable{:,columnhpred}=hpred;
            else
                fTable=[fTable table(hpred)];
            end
            PpredX=Penc(:,1);PpredY=Penc(:,2);
            columnPpredX=find(ismember(fieldnames(fTable),'PpredX'));
            if columnPpredX
                fTable{:,columnPpredX}= PpredX;
            else
                fTable=[fTable table( PpredX)];
            end
            columnPpredY=find(ismember(fieldnames(fTable),'PpredY'));
            if columnPpredY
                fTable{:,columnPpredY}= PpredY;
            else
                fTable=[fTable table( PpredY)];
            end

            %fTable1= [fTable table(hpred,PpredX,PpredY)];
            writetable(fTable,fileName,'Sheet',sheetName);

        end

        function [Te,Pe, EncodingDist_TP,EncodingAng_alpha,he]= predict_P(self,a,b, c,d,O,OT, TP, v_OT, alpha)
             [Te,Pe, EncodingDist_TP,EncodingAng_alpha,he]=self.encodingErrorLocation(a,b, c,d, O,OT, TP, v_OT, alpha);
        end 
        
        function Erem_O = predict_O(self, p,q,m,n,Obj,P,h,Ppred,hpred)
             %[Te,Pe, EncodingDist_TP,EncodingAng_alpha]=encodingError1(a,b, c,d, O,OT, TP, v_OT, alpha);

             %VectorCal=collectiveVectorCal;
%              T= v_OT + O;
%              v_TP=P-T;
%              h=self.VectorCal.getangle_vAtovB([0,1],v_TP);

             Predict_v1=self.responseErrorLocation(p,q,m,n, hpred, Ppred, Obj);
             %Predict_v1=self.responseError1(p,q,m,n, hpred, Ppred, Obj);
             %VectorCal=collectiveVectorCal;
             Predict_v=self.VectorCal.RotateVecArrayLeft(h-hpred,Predict_v1);%he- h???
             Erem_O =  Predict_v + P;  %%% predicted pointing O'
        end 
        
        function simu_Oresp=getSiOresp(self,P,O,alpha,OT,TP,v_OT,simulationPara,sheetName,WhichObjects,fileName)
            %simulationPara has 9 parameters, the last one is the noise
            %mse for x and y, which should be mse of location residuals/2
            
            a=simulationPara(1);
            b=simulationPara(2);
            c=simulationPara(3);
            d=simulationPara(4);
            p=simulationPara(5);
            q=simulationPara(6);
            m=simulationPara(7);
            n=simulationPara(8);
            noiseVar=simulationPara(9); % which is the variance of x or y, one dimension, so two dimension should be 10
%             sum1=a+b+c+d;
%             sum2=p+q+m+n;
%             if sum1==0
%                 a=1;b=0;c=1;d=0;
%             end


            for i=1:numel(O(1,:))/2
                 Erem_O(:,2*i-1:2*i)=repmat(NaN,numel(O(:,1)),2);
                 index=find(~isnan(O(:,2*i-1)));
                 Erem_O(index,2*i-1:2*i) = self.predict_Location(a,b, c,d, p,q,m,n,O(index,1:2),O(index,2*i-1:2*i),P(index,:), OT(index), TP(index), v_OT(index,:), alpha(index));
                 %error =[error;(Oresp(index,2*i-1)-Erem_O(index,1)).^2+(Oresp(index,2*i)-Erem_O(index,2)).^2];
                 %index=find(isnan(O(:,2*i-1)));    
            end
            simu_Oresp=Erem_O;
            
            %on1/20/2022 I will add some noise to the simu_Oresp
            if isempty(self.noise)
            %if ~ismember(sheetName,keys(self.noise))
                %self.noise(sheetName)=randn([numel(simu_Oresp(:,1)) numel(simu_Oresp(1,:))])*sqrt(noiseVar);
                self.noise=randn([numel(simu_Oresp(:,1))*2 (self.objectTotal+1)*2])*sqrt(noiseVar);
            end
            
            fTable= readtable(fileName,'Sheet',sheetName);
            ID=fTable.ID-640;
            noise1=self.noise;
            noiseAll=self.noise(ID,:);
            
            noiseSelectedAll=[];
            for i=1:length(WhichObjects)
                column=2*WhichObjects(i)+1;
                noiseSelectedAll=[noiseSelectedAll,noiseAll(:,column:column+1)];          
            end
            
            simu_Oresp=simu_Oresp+noiseSelectedAll;
        end
        
        function [PprimeBI, hprimeBI]=getPprimehprimeBi(self,fileName,sheetName,objectTotal)  
            
            whichObjects=[0:objectTotal];
            [P,T,O,Oprime,alpha,OT,TP,v_OT,Pprime,hprime,h,Ppred,hpred]=self.getPathInfo(fileName,sheetName,'all',self.objectTotal,whichObjects,self.legRatio);%all O and X
            
            if ~isempty(self.simulationPara) %do simulation instead
                Oprime= self.getSiOresp(P,O,alpha,OT,TP,v_OT,self.simulationPara,sheetName, whichObjects,fileName);
            end
            
            %vectorCal=collectiveVectorCal;
%             v_TP=P-T;
%             h=vectorCal.getangle_vAtovB([0,1],v_TP);
            if objectTotal==2
                LastRowofFiveObjects=0;
            else
                LastRowofFiveObjects=640;
            end
            for i=1:LastRowofFiveObjects% numel()
                X=[];
                Xprime=[];
                for m=1:objectTotal
                    X=[X;O(i,2+2*m-1:2+2*m)];
                    Xprime=[Xprime;Oprime(i,2+2*m-1:2+2*m)];
                    %X=[O(i,3:4);O(i,5:6);O(i,7:8);O(i,9:10)];
                end
                    
                
                %Xprime=[Oprime(i,3:4);Oprime(i,5:6);Oprime(i,7:8);Oprime(i,9:10)];
                output=self.HprimePrimeForOnePathFunction_bi_regression (O(i,1:2),X, P(i,:), h(i), Oprime(i,1:2), Xprime);
                hprime(i,1)=output.hprime; %column BC
                Pprime(i,:)=output.Pprime;
                Rsq(i,:)=output.Rsq;
            end

            for i=LastRowofFiveObjects+1:numel(h)
                X=[O(i,3:4);O(i,5:6)];
                Xprime=[Oprime(i,3:4);Oprime(i,5:6)];
                output=self.HprimePrimeForOnePathFunction_bi_regression (O(i,1:2),X, P(i,:), h(i), Oprime(i,1:2), Xprime);
                hprime(i,1)=output.hprime; %column BC
                Pprime(i,:)=output.Pprime;
                Rsq(i,:)=output.Rsq;
            end

            %fTable= readtable(fileName,'Sheet',sheetName);
            hprimeBI=hprime;
            PprimeBI=Pprime;
            
            
        end
        
        function rec=HprimePrimeForOnePathFunction_bi_regression (self,O, XY, P, h, Oprime, XYprime)

            myData=[[O;XY], [Oprime;XYprime]];
            

            output=bi_regressV3(myData(:,3),myData(:,4),myData(:,1),myData(:,2),P(:,1),P(:,2)); %AB as independent
            hprime=h-output.regression(5); %column BC
            %O=[myData1(1,1), myData1(1,2)]
            Pprime=output.dependent;
        

            rec.Pprime=Pprime; rec.hprime=hprime;
            rec.Rsq=output.regression(7);
       
        end
        
        function [mdl_ang,p_ang]=scatterPlotAngle(self,Ax,Ay)%Ax will be changed to positive, Ay will change the sign if Ax changes the sign.
            %mdl_dis=self.regressionLine(v_PprimeO_dAll, v_PprimeOresp_dAll);
                   
                   
                    
                    Ay=Ay.*Ax./abs(Ax);
                    Ax=abs(Ax);

                    Ay=self.VectorCal.normalizeAngBMinimizeDisToA(Ax,Ay);%VectorCal.ChangeAngleTo_0to360(Angle_TPOresp1);
                    
                    mdl_ang=self.regressionLine(Ax,Ay);
                    p_ang=polyfit(Ax,Ay,1);
        end
        
        
        function visualizePointsPairs(self,P1,P2)
            
            x1=P1(:,1);
            y1=P1(:,2);
            x2=P2(:,1);
            y2=P2(:,2);

        %     col=zeros(20,3)
            %p=colormap
        %     for i=1:20
        %       a=y(i);
        %       col(i,:)=p(a,:);
        %       stem3(i,i,x(i),'Color',col(i,:))
        %       hold on
        %     end
%             AxUni=unique(P1);
%             ColorUni=[1:round(256/numel(AxUni)):(numel(AxUni)-1)*round(256/numel(AxUni))+1]
%             ColorUni=ColorUni(randperm(length(ColorUni)))

            for i=1:numel(x1)
                 col=zeros(1,3);
                 p=colormap;
                 %a=find(AxUni==Angle1(i));
%                 %a=round(Angle1(i)/max(Angle1)*256);
                col=[p(1,:);p(2,:)];
                %plot([x1(i),x2(i)],[y1(i),y2(i)]);               
                plot(x1(i),y1(i),'o');
                plot(x2(i),y2(i),'x');
                plot([x1(i),x2(i)],[y1(i),y2(i)]);
                hold on
            end

        end
        
        function  [O,Oprime,PprimeBI, hprimeBI]=getPpredhpredFlying(self,fileName,sheetName)
            [PprimeBI, hprimeBI]=self.getPprimehprimeBi(fileName,sheetName,self.objectTotal); 
        end
        
        function [Ppred,hpred]=getPpredhpred(self,sheetName,O, OT, TP, v_OT, alpha)
            [P,T,O,Oresp,alpha,OT,TP,v_OT,Pprime,hprime,h,Ppred,hpred]=self.getPathInfo(fileName,sheetName,flag,self.objectTotal,self.whichObjects,self.legRatio);
            
            %self.encoParaTable=unique(self.encoParaTable);
            %self.encoParaTable.RowNames = self.encoParaTable.sheetName;
            
            ind=find(self.encoParaTable.sheetName==sheetName);
            para=self.encoParaTable.encoPara(ind(1),:);
            
            [Tenc,Ppred, EncodingDist_TP,EncodingAng_alpha,hpred]= self.predict_P(para(1),para(2), para(3),para(4), O, OT, TP, v_OT, alpha);
            %PpredX=Penc(:,1);PpredY=Penc(:,2);    
        end 
        
        function self=setobjectTotal(self,x)
            self.objectTotal=x;
            %self.IterationNumber=5;
        end
        
        function [dx,dy,Ax,Ay]=dis_angPredResphomingV(self,O,P,h,Oresp,Pprime,hprime)
                    v_PprimeO_dAll=[];
                    v_POresp_dAll=[];
                    Angle_hprimePprimeOAll=[];
                    Angle_hPOrespAll=[];
                    for i=1:numel(O(1,:))/2
                        index=find(~isnan(O(:,2*i-1)));

                        v_PprimeO=self.calculate_vector(Pprime(index,:),O(index,2*i-1:2*i));
                        [v_PprimeO_d, tem]=self.VectorCal.getBearingofAB(Pprime(index,:),O(index,2*i-1:2*i));

                        %v_TPprime=calculate_vector(T(index,:),Pprime(index,:));
                        %v_OT=calculate_vector(O,T);
                        v_POresp=self.calculate_vector(P(index,:),Oresp(index,2*i-1:2*i));
                        
                        [v_POresp_d, tem]=self.VectorCal.getBearingofAB(P(index,:),Oresp(index,2*i-1:2*i));
                        v_PprimeO_dAll=[v_PprimeO_dAll;v_PprimeO_d];
                        v_POresp_dAll=[v_POresp_dAll;v_POresp_d];

                        v_hprime=[sind(hprime(index)),cosd(hprime(index))];
                        v_h=[sind(h(index)),cosd(h(index))];

                        Angle_hprimePprimeO=self.VectorCal.getangle_vAtovB(v_hprime,v_PprimeO);
                        Angle_hPOresp=self.VectorCal.getangle_vAtovB(v_h,v_POresp);
                        Angle_hprimePprimeOAll=[Angle_hprimePprimeOAll;Angle_hprimePprimeO];
                        Angle_hPOrespAll=[Angle_hPOrespAll;Angle_hPOresp];

                    end
            
                   
                    dx=v_PprimeO_dAll;
                    dy=v_POresp_dAll;
                    Ax=Angle_hprimePprimeOAll;
                    Ay=Angle_hPOrespAll;
                   
                    
                    Ay=Ay.*Ax./abs(Ax);
                    Ax=abs(Ax);

                    % figure
                    % plot(Angle_TPprimeO1All, Angle_TPprimeOresp1All,'.');
                    % axis equal
                    Ay=self.VectorCal.normalizeAngBMinimizeDisToA(Ax,Ay);%VectorCal.ChangeAngleTo_0to360(Angle_TPOresp1);
        end 
            
        function [OAll, OrespAll]=selectOOresp(self,O,Oresp,WhichObjects)%WhichObjects=[0,1,3] will get O, X1, X3
                OAll=[];
                OrespAll=[];
                for i=1:length(WhichObjects)
                    column=2*WhichObjects(i)+1;
                    OAll=[OAll,O(:,column:column+1)];
                    OrespAll=[OrespAll,Oresp(:,column:column+1)];
                end
        end
        
        
        function [OAll,OrespAll,Erem_OAll]=getORPAllOriginalFormat(self,x,O,P,OT,TP,Oresp,v_OT,alpha,Pprime,hprime,h,Ppred,hpred,whichModel,flag) % do not append below vertically
            
              %whichModel=1, encoding only, 2, response only, 3 combined
             if whichModel==1
                 a=x(1);
                 b=x(2);
                 c=x(3);
                 d=x(4);
                 p=1;
                 q=0;
                 m=1;
                 n=0;
             elseif whichModel==2
                 a=1;
                 b=0;
                 c=1;
                 d=0;
                 p=x(1);
                 q=x(2);
                 m=x(3);
                 n=x(4);
             else
                 a=x(1);
                 b=x(2);
                 c=x(3);
                 d=x(4);
                 p=x(5);
                 q=x(6);
                 m=x(7);
                 n=x(8);
             end

              OrespAll=[];
              Erem_OAll=[];
              OAll=[];

              %modelLogL=0;
              %baseLogL=0;
             for i=1:numel(O(1,:))/2
                 index=find(~isnan(O(:,2*i-1)));
                 if strcmp(flag,'all')  %step 2 always uses 'all'
                    
                     Erem_O = self.predict_Location(a,b, c,d, p,q,m,n,O(index,1:2),O(index,2*i-1:2*i),P(index,:), OT(index), TP(index), v_OT(index,:), alpha(index));

                     Erem_OAll=[Erem_OAll,Erem_O];
                     OrespAll=[OrespAll,[Oresp(index,2*i-1),Oresp(index,2*i)]];
                     OAll=[OAll,[O(index,2*i-1),O(index,2*i)]];                 
                 elseif strcmp(flag,'enco') %only for step 3
                    [Tenc,Penc, EncodingDist_TP,EncodingAng_alpha,henc] = self.predict_P(a,b, c,d, O(index,1:2), OT(index), TP(index), v_OT(index,:), alpha(index));
                    Erem_O =Penc;
                    Erem_OAll=[Erem_OAll,Erem_O];
                    OrespAll=[OrespAll,[Pprime(index,2*i-1),Pprime(index,2*i)]];
                    OAll=[OAll,[P(index,1),P(index,2)]]; 
                    %error =[error;(Pprime(index,2*i-1)-Ppred(index,1)).^2+(Pprime(index,2*i)-Ppred(index,2)).^2];
                 else %only for step 3
                     Erem_O = self.predict_O(p,q,m,n,O(index,2*i-1:2*i),P(index,:), h(index),Ppred(index,:),hpred(index));
                     Erem_OAll=[Erem_OAll,Erem_O];
                     OrespAll=[OrespAll,[Oresp(index,2*i-1),Oresp(index,2*i)]];
                     %base o should be affected by pe he already
                     Base_O=self.predict_O(1,0,1,0,O(index,2*i-1:2*i),P(index,:), h(index),Ppred(index,:),hpred(index));                 %OAll=[OAll;[O(index,2*i-1),O(index,2*i)]]; 
                     OAll=[OAll,Base_O];
                 end

             end
         end  
         function visualizeAngles(self,r,Angle,c)


            fg=figure;
            r1=r;
            %r2=2;
            self.drawCircle(fg,[0,0], r1);

            for k=1:length(Angle(1,:))
                Angle1=self.VectorCal.RightToLeftAngle(Angle(:,k)); %%Angle2=self.VectorCal.RightToLeftAngle(Angle2);
                x1=r1*cosd(Angle1);
                y1=r1*sind(Angle1);

                                  
                plot(x1,y1,c{k});
                hold on                 
            end

         end
        
        function drawCI(self,r0,cM,cSD,c)
            center=[0,0];
            low=(cM-cSD);
            high=(cM+cSD);
            for i=1:length(cM)
                r=r0(i);
                theta1 = linspace(low(i),high(i));
                theta=self.VectorCal.RightToLeftAngle(theta1);
                xc = r*cosd(theta)+center(1);
                yc = r*sind(theta)+center(2);
                plot(xc,yc,c{i});
                hold on
                axis equal
            end
        end 
        
        function [t,p]=Dietterich(self,dabAll)%2by5, 10 d in a column vector
    
            vardabAll=self.getvardabAll(dabAll);
            t=dabAll(1,:)./sqrt(sum(vardabAll)/5);

            p=1-tcdf(abs(t),5);

        end
        function [F,p]=Alpaydin(self,dabAll)
            vardabAll=self.getvardabAll(dabAll);
            F=sum(dabAll.^2)./(2*sum(vardabAll));
            p=1-fcdf(F,10,5);
        end
        function vardabAll=getvardabAll(self,dabAll)
            vardabAll=[];
            for i=1:length(dabAll(:,1))/2
                a=dabAll(2*i-1,:);
                b=dabAll(2*i,:);
                meanab=mean([a;b]);
                vardab=(a-meanab).^2+(b-meanab).^2;
                vardabAll=[vardabAll;vardab];
            end
        end
        
        function lr=pairLr(self,logL)
            dlogLAll=[];
            for i=1:length(logL)
                for j=i+1:length(logL)
                    dlogL=logL(i)-logL(j);
                    dlogLAll=[dlogLAll,dlogL];
                end     
            end
            lr=exp(dlogLAll);
        end
        
        
        function pRsqAll=pairpRsq(self,mseMat)
             pRsqAll=[];
            for i=1:length(mseMat(1,:))
                for j=i+1:length(mseMat(1,:))
                    pRsq=1-mseMat(:,j)./mseMat(:,i);
                    pRsqAll=[pRsqAll,pRsq];
                end     
            end
            
        end
        
        
        function [d,A]=dis_ang_homingV(self,O,P,h) %homing length and angle
                    v_PO_dAll=[];
         
                    Angle_hPOAll=[];
                    
                    for i=1:numel(O(1,:))/2
                        index=find(~isnan(O(:,2*i-1)));

                        v_PO=self.calculate_vector(P(index,:),O(index,2*i-1:2*i));
                        [v_PO_d, tem]=self.VectorCal.getBearingofAB(P(index,:),O(index,2*i-1:2*i));

                        
                        v_PO_dAll=[v_PO_dAll;v_PO_d];
                        
                        v_h=[sind(h(index)),cosd(h(index))];

                        Angle_hPO=self.VectorCal.getangle_vAtovB(v_h,v_PO);
                        
                        Angle_hPOAll=[Angle_hPOAll;Angle_hPO];
                       
                    end
            
                   
                    d=v_PO_dAll;
                    
                    A= Angle_hPOAll;
                   
        end 
        
    end
    
end