classdef proj39collection
    properties (SetAccess = private)
        VectorCal=collectiveVectorCal;
    end
    methods 
        function self=proj39collection()  
            %newpath='\\weiminmlab.synology.me\labmeeting\labmeeting\circular statistics\circular statistic for matlab\CircStat2012a'
            newpath='\\129.128.52.128\labmeeting\labmeeting\circular statistics\circular statistic for matlab\CircStat2012a'
            userpath(newpath);
        end
        function [AxUni, Aymean]=NamesForUnqueX(self,Ax,Ay)
            AxUni=unique(Ax);
            Aymean=cell(size(AxUni,1),1);
            for i=1:numel(AxUni)
                Ay1=Ay(find(Ax==AxUni(i)));%strcmp(fTable.env_building,Uni(i))
                Aymean1=Ay1{1,:};%mean(Ay1);
                Aymean{i}=Aymean1;
            end
            %AyCmean
        end
        function [AxUni, Aymean]=ValueForUnqueX(self,Ax,Ay)
            AxUni=unique(Ax);
            Aymean=[];%cell(size(AxUni,1),1);
            for i=1:numel(AxUni)
                Ay1=Ay(find(strcmp(Ax,AxUni{i,:})));
                Aymean1=Ay1(1);%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 [meanOresp,unicateM]=getMeanORP1(self,Oresp,cateM)
         %[meanM,unicateM]=getMeanORP([Oresp,alpha,OT,TP,v_OT,Pprime,hprime,h,Ppred,hpred],[expNo,P,T,O]);
                A=Oresp;
                
                [C,ia,ic] = unique(cateM);
                B=[];
                for i=1:length(ia) %((ia) stores unique categoray)
                    index=ic==i;%get category i
                    B=[B;mean(Oresp(index,:))];       
                end
                meanOresp=B(:,1:length(Oresp(1,:)));
                unicateM=C;
        
                
        %         O=B(:,1:2);
        %         Oresp=B(:,3:4);
        %         Opred1=B(:,5:6);
        %         Opred2=B(:,7:8);
        %         Opred3=B(:,9:10);
        end
        function [meanOresp,unicateM]=getCircularMeanORP1(self,Oresp,cateM)
         %[meanM,unicateM]=getMeanORP([Oresp,alpha,OT,TP,v_OT,Pprime,hprime,h,Ppred,hpred],[expNo,P,T,O]);
                A=Oresp;
                
                [C,ia,ic] = unique(cateM);
                B=[];
                for i=1:length(ia) %((ia) stores unique categoray)
                    index=ic==i;%get category i
                    
                    B=[B;circMean(Oresp(index,:))];       
                end
                meanOresp=B(:,1:length(Oresp(1,:)));
                unicateM=C;
        
                
        %         O=B(:,1:2);
        %         Oresp=B(:,3:4);
        %         Opred1=B(:,5:6);
        %         Opred2=B(:,7:8);
        %         Opred3=B(:,9:10);
        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);
                   
                   VectorCal=collectiveVectorCal;
                    
                    Ay=Ay.*Ax./abs(Ax);
                    Ax=abs(Ax);
        
                    Ay=VectorCal.normalizeAngBMinimizeDisToA(Ax,Ay);%VectorCal.ChangeAngleTo_0to360(Angle_TPOresp1);
                    
                    mdl_ang=self.regressionLine(Ax,Ay);
                    p_ang=polyfit(Ax,Ay,1);
        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,180]);           %xlim([-10,30]); xlim([-3,59]);
                ylim([-20,180]);           %ylim([-30,60]); ylim([-30,80])
                xlabel('Predicted angle errors') 
                ylabel('Mean response angle erros') 
            
                 mdl=mdl_angle;
        end
        
        function [mdregression1,tb1]=myregression (self,T, VarName)
            tb1=table()
            for i=1:size(VarName,2)
                tb1.(VarName{i})=T.(VarName{i});
            end
            mdregression1 = fitlm(tb1,'ResponseVar',VarName{1});
        end 
        
        function d=earthD(self,A,B) %A,B in degree
            %https://www.omnicalculator.com/other/latitude-longitude-distance
            % θ₁, φ₁ – First point latitude and longitude coordinates;
            % θ₂, φ₂ – Second point latitude and longitude coordinates;
            % R – Earth's radius (R = 6371 km); and
            % d – Distance between them along Earth's surface.
            %d = 2R × sin⁻¹(√[sin²((θ₂ - θ₁)/2) + cosθ₁ × cosθ₂ × sin²((φ₂ - φ₁)/2)])
            R = 6371000 %m
            A=A*pi()/180;
            B=B*pi()/180;
            theta1=A(:,1);  phi1=A(:,2);
            theta2=B(:,1);  phi2=B(:,2);
           
            d= 2*R*asin( sqrt( (sin((theta2-theta1)/2))^2 + cos(theta1)*cos(theta2)*(sin((phi2-phi1)/2))^2 ) );
        end
        
        
        function theta = earthbearing(self,A,B)
              
                lat1=A(:,1);
                long1=A(:,2);
                lat2=B(:,1);
                long2=B(:,2);
        
                %   % MATLAB
                long1 = long1 * (pi/180);
                long2 = long2 * (pi/180);
                lat1 = lat1 * (pi/180);
                lat2 = lat2 * (pi/180);
                
                theta = atan2(sin(long2 - long1) * cos(lat2), cos(lat1) * sin(lat2) - sin(lat1) * cos(lat2) * cos(long2 - long1));
                theta = theta * 180 / pi;
                
                if theta < 0
                    theta = theta + 360;
                end
        
        end
        
        function [x,y,d,theta,d1]=XYfromLatLong(self,referenceLocation,B)
            A=referenceLocation;
            theta = proj39collection.earthbearing(A,B);
            d=proj39collection.earthD(A,B);
            [x,y]=proj39collection.getxyfromdisandbearing(d,theta);
            d1=sqrt(x^2+y^2);    
            %get x
            
        
            % A=referenceLocation;
            % BxAy=[B(:,1),A(:,2)];
            % AxBy=[A(:,1),B(:,2)];
            % x=earthD(A,BxAy);
            % y=earthD(A,AxBy);
            % d=sqrt(x^2+y^2)
        
        
            %get y
        
        end
        function [x,y]=getxyfromdisandbearing(self,dis,bearing)%bearing in degree
            %bearing in terms of z positive or north
            x=dis*sind(bearing);
            y=dis*cosd(bearing);
        end
        
        function T=pilotTable(self,fTable, unicateM, varNames,method)
            tbCate=table();
            for i=1:size(unicateM,2)
                tbCate.(unicateM{i})=fTable.(unicateM{i});
            end
        
            % tbCate.indoor_outdoor=fTable.indoor_outdoor;
            % %tbCate.Enviornmentname=fTable.Enviornmentname;
            % tbCate.participant_number=fTable.participant_number;
            A=[];
            for i=1:size(varNames,2)
                A=[A,fTable.(varNames{i})];
            end
            if strcmp(method,'mean')
                [meanM,unicateM]=self.getMeanORP1(A,tbCate);
            elseif strcmp(method,'circMean')
                [meanM,unicateM]=self.getCircularMeanORP1(A,tbCate);
            else
                print "wrong method"
            end 
            
            %[meanM,unicateM]=getMeanORP1([fTable.abs_pointerror,fTable.PTT,fTable.MRT,fTable.continuous_unfam_fam,fTable.euclidian_distance,fTable.campus_experience,fTable.absheadingError,fTable.distanceError,fTable.leastErrSquared],tbCate);
            %[meanM,unicateM]=getMeanORP([fTable.abs_pointerror,fTable.PTT,fTable.MRT,fTable.euclidian_distance,fTable.campus_experience,fTable.envrionmentFamilarity],[fTable.indoor_outdoor,fTable.groupenvrionmentFamilarity,fTable.participant_number]);
            %A=[meanM,unicateM];
            
            T = array2table(meanM,...
                'VariableNames',varNames);
            
            T=[unicateM,T];
        end 
        
        function TA=reshapetoRepeated(self,Ts,varNames,levels,newvarNames)
            TA=Ts(1:size(Ts(:,1),1)/levels, :);
            
            for i=1:size(varNames,2)
                varname=varNames{i};
                A=Ts.(varname);
                B=reshape(A,size(A,1)/levels,levels);
                tbB=array2table(B, 'VariableNames',newvarNames{i});
                TA=[TA, tbB];
            end
            
        end
        function output=circMean(self,A)
        
            output= circ_rad2ang(circ_mean(circ_ang2rad(A)));
        
        end

        function output=circConfmean(self,A)
        
            output= circ_rad2ang(circ_confmean(circ_ang2rad(A)));
        
        end
       
        
        function output=bidiRegressWithoutTransformation(predictor, predicted)
            meanpredicted=mean(predicted);
            meanpredicted = repmat(meanpredicted,size(predicted(:,1),1));
            distanceUnPredicted=distance(predictor,predicted);
            distanceTotal=distance(meanpredicted,predicted);
            output.meandistanceUnPredicted=mean(distanceUnPredicted);
            output.meandistanceTotal=mean(distanceTotal);
            output.explained=1-distanceUnPredicted/distanceTotal;
        end
        
        function  d=distance(self,predictor,predicted)
            d=sqrt(sum((predictor-predicted).^2,2));
        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 visualOneBuilding(P1,P2)
        %     % %  fields=fieldnames(P);
        %     % % c={'g','r','b','k'};
        %     % % figure
        %     % % for i=1:length(fields)
        %     % %     XY0=P.(fields{i});
        %     % %     if scatterFlag
        %     % %         scatter(XY0(:,1),XY0(:,2),['.',c{i}]);
        %     % %         scale=0.95;
        %     % %     else
        %     % %         scale=2;
        %     % %     end
        %     % %     hold on
        %     % %     if residualFlag
        %     % %         output.(fields{i})=ellipsePandE_Residual20210716(1, scale,XY0,['+',c{i}]);
        %     % %     else
        %     % %         output.(fields{i})=ellipsePandE20210716(1, scale,XY0,['+',c{i}]);
        %     % %     end
        %     % %     hold on
        %     % end
        %     %output.(fields{i})=ellipsePandE_Residual20210716(1, scale,XY0,['+',c{i}]);
        %     scale=0.95;
        %     XY0=P2.(fields{1})
        %     ellipsePandE_Residual20210716(1, scale,XY0,['+',c{i}]);
        % 
        %     plot(P1.(fields{1}(1)),P1.(fields{1}(2)),'o')
        % 
        %     title(sheetName);
        %     axislim=10;
        %     xlim([-axislim axislim])
        %     ylim([-axislim axislim])  
        %     grid on
        %     pbaspect([1 1 1])
        % 
        % end
        function angle = convertangleto360(self,x)
                angle = mod(x,360);
                % if x < 0
                %     angle = x +360;
                % elseif x > 360
                %     angle = x-360;
                % else
                %     angle = x;
        end

        function bearing = getbearingbetween(self,v1,v2)
            bearing = acos(dot(v1,v2)/(norm(v1)*norm(v2)))*180/pi;
            crossproduct = cross(v1,v2);
            dotproduct = dot(crossproduct, [0,1,0]);
            if dotproduct > 0
                bearing = -bearing +360;
            end
        end
        function bearing = getbearingA (self,x,A,i)

             P1A = [A(1) - x(1),0, A(3) - x(2)];
             P1Zplus = [0,0,1]; 
            
             bearingA_raw = getbearingbetween(P1A,P1Zplus);
             bearing = convertangleto360(bearingA_raw - i);
        end
        function angle = convertangleto180(self,x)
            x = mod(x,360); %change to 0 to 360
            index=x>180;
            x(index)=x(index)-360;
            angle=x;
            % if x < -180
            %     angle = x +360;
            % elseif x > 180
            %     angle = x-360;
            % else
            %     angle = x;
        end


        function [mse,output,logL]=visualizeORP(proj36,outputfileName,fileName,sheetName,half,scatterFlag,residualFlag,meanFlag)
            
        
            [O,Oresp,Opred1]=getORP(proj36,outputfileName,'x1',fileName,sheetName,half,1,meanFlag); 
            [O,Oresp,Opred2]=getORP(proj36,outputfileName,'x2',fileName,sheetName,half,2,meanFlag);
            [O,Oresp,Opred3]=getORP(proj36,outputfileName,'x3',fileName,sheetName,half,3,meanFlag);
            
        %     if meanFlag
        %         %[O,Oresp,Opred1,Opred2,Opred3]=getMeanORP(O,Oresp,Opred1,Opred2,Opred3);
        %         [meanOresp,unicateM]=getMeanORP(Oresp,[O,Opred1,Opred2,Opred3]);
        %         Oresp=meanOresp;
        %         O=unicateM(:,1:2);
        %         Opred1=unicateM(:,3:4);
        %         Opred2=unicateM(:,5:6);
        %         Opred3=unicateM(:,7:8);
        %     end
            
            
            P.XY0=O-Oresp;
            mse.xy0=getXYMSE(P.XY0);
            logL.xy0=proj36.MaxLikelihoodLocation(O,Oresp);
            
            P.XY1=Opred1-Oresp;  
            mse.xy1 =getXYMSE(P.XY1);
            logL.xy1=proj36.MaxLikelihoodLocation(Opred1,Oresp);
            
            P.XY2=Opred2-Oresp;
            mse.xy2 =getXYMSE(P.XY2);
            logL.xy2=proj36.MaxLikelihoodLocation(Opred2,Oresp);
            
            P.XY3=Opred3-Oresp;
            mse.xy3 =getXYMSE(P.XY3);
            logL.xy3=proj36.MaxLikelihoodLocation(Opred3,Oresp);
            
            
            fields=fieldnames(P);
            c={'g','r','b','k'};
            figure
            for i=1:length(fields)
                XY0=P.(fields{i});
                if scatterFlag
                    scatter(XY0(:,1),XY0(:,2),['.',c{i}]);
                    scale=0.95;
                else
                    scale=2;
                end
                hold on
                if residualFlag
                    output.(fields{i})=ellipsePandE_Residual20210716(1, scale,XY0,['+',c{i}]);
                else
                    output.(fields{i})=ellipsePandE20210716(1, scale,XY0,['+',c{i}]);
                end
                hold on
            end
            
            plot(0,0,'o')
            
            title(sheetName);
            axislim=10;
            xlim([-axislim axislim])
            ylim([-axislim axislim])  
            grid on
            pbaspect([1 1 1])
        end 

        function visualizeAnglePairs(self,Angle1,Angle2)
           
            %Angle1=self.VectorCal.RightToLeftAngle(Angle1); Angle2=self.VectorCal.RightToLeftAngle(Angle2);
            fg=figure
            r1=1
            r2=1.2
            fg=visualizeAngle(self,Angle1,r1,[0 0],fg,'r')
            fg=visualizeAngle(self,Angle2,r2,[0 0],fg,'b')

            text(0,0+r2+0.3,'0')
            text(0+r2+0.1,0,'90')
            text(0,0-r2-0.1,'180')
            text(0-r2-0.3,0,'-90')
            xlim([-2 2])
            ylim([-2 2])
            set(gca,'XTick',[],'YTick',[])
            
        end
        function fg=visualizeAngle(self,Angle,r,location,fginput, color)
           
            ci1=self.circConfmean(Angle);
          
            %ci2=self.circConfmean(Angle2);

            mean=self.circMean(Angle);
            meanl=mean-ci1;
            meanr=mean+ci1;
            
            sign=meanl*meanr;
            if sign>0
                symble='*'
            else
                symble='ns'
            end

            mean=self.VectorCal.RightToLeftAngle(mean);
            meanl=self.VectorCal.RightToLeftAngle(meanl);
            meanr=self.VectorCal.RightToLeftAngle(meanr);

            meanx=r*cosd(mean);
            meany=r*sind(mean);
           
            Angle1=self.VectorCal.RightToLeftAngle(Angle); %Angle2=self.VectorCal.RightToLeftAngle(Angle2);
            fg=figure(fginput)
            plot([location(1),location(1)+meanx],[location(2),location(2)+meany],color);
            hold on
            self.drawArc(fg,location, r+0.1, meanl, meanr,color);
            text(location(1)+(1+0.1)*meanx,location(2)+(1+0.1)*meany,symble)
            hold on
            %r1=r;
            %r2=2;
            self.drawCircle(fg,location, r,color);
            %self.drawCircle(fg,[0,0], r2);
            x1=r*cosd(Angle1);
            y1=r*sind(Angle1);
            %x2=r2*cosd(Angle2);
            %y2=r2*sind(Angle2);

            scatter(x1,y1);
            hold on

        end
        function drawCircle(self,fg,center, r,color)
            theta = linspace(-pi,pi);
            xc = r*cos(theta)+center(1);
            yc = r*sin(theta)+center(2);
            plot(xc,yc,color);
            hold on
            axis equal
        end 
        function drawArc(self,fg,center, r,leftAngle,rightAngle, color)
            theta = linspace(leftAngle/180*pi,rightAngle/180*pi);
            xc = r*cos(theta)+center(1);
            yc = r*sin(theta)+center(2);
            plot(xc,yc,color);
            hold on
            axis equal
        end 

         function [FormStr, mdl_angle]=regressionLineWithoutFigure(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);
                % x=2.5
                % y=mdl_angle.Coefficients.Estimate(2).*x+mdl_angle.Coefficients.Estimate(1)
                % text(x,y, FormStr)
                
                %legend({'','Diagonal line',FormStr, 'Location','southeast','FontSize','14'})
                %legend({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,180]);           %xlim([-10,30]); xlim([-3,59]);
            %     ylim([-20,180]);           %ylim([-30,60]); ylim([-30,80])
            %     xlabel('Predicted angle errors') 
            %     ylabel('Mean response angle erros') 
            
                % mdl=mdl_angle;
        end
    end
end