
%
clear all
close all

HotelingT=collectiveHotelingTsquared

T12PositionInfo= readtable('T12PositionInfo20240113.csv');
unicateM=table(T12PositionInfo.indoor_outdoor,T12PositionInfo.Enviornmentname);
[C,ia,ic] = unique(unicateM);
%B=[];
%%compare prime with correct, testing bias
tprime2c=table();
for i=1:length(ia) %((ia) stores unique categoray)
    index=ic==i;%get category i
    key=[C.Var2{i}(1:2),num2str(C.Var1(i))];
    key = erase(key," ");
    %key=key(1:5);
    Pcorrect.(key)=[T12PositionInfo(index,:).('positionX'), T12PositionInfo(index,:).('positionY')];
    Pprime.(key)=[T12PositionInfo(index,:).('positionprimeX'), T12PositionInfo(index,:).('positionprimeY')];

    %P.{[]}
    %B.=[B;mean(Oresp(index,:))];       
end

figure 
fields=fieldnames(Pprime);
c=jet(length(fields));
% col=zeros(1,3);
cp=colormap;
%c=p(int16(256/length(fields))*i+1,:);
%r = randi([1 256],12,1)
for i=1:length(fields)  
    figure
    c=cp(int16(256/length(fields))*i+1,:)
    %p(int16(256/length(fields))*i+1,:);
     %a=find(AxUni==Angle1(i));
    %                 %a=round(Angle1(i)/max(Angle1)*256);
    %c=[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)]);
    %colororder(c);
   
    %c={'g','r','b','k'};
    thiscolor=c;
    P1=Pcorrect.(fields{i});
    P2=Pprime.(fields{i});
    %scatter(P2(:,1),P2(:,2),['.',c{1}]);
    %scatter(P2(:,1),P2(:,2),'.','ColorVariable',c{1})
    %plot(P2(:,1),P2(:,2),'.','Color',c{1})
    plot(P2(:,1),P2(:,2),'.','Color',thiscolor)
    hold on
    %ellipsePandE20210716(1, 0.95,P2,['+',c{1}])
    output.(fields{i})=ellipsePandE20231125(2, 0.95,P2,'+',thiscolor);
    hold on
    plot(P1(1,1),P1(1,2),'o','Color',thiscolor);
    hold on
    text(P1(1,1),P1(1,2),fields{i},'Color','b');
    Mu=output.(fields{i}).Mu;
    text(Mu(1),Mu(1,2),[fields{i},char(39)],'Color','k');
    %plot(P1(1,1),P1(1,2),['o',c{1}]);
    
    %do hythothesis test between prime and correct positions using
    %HotelingT
    M=[P2,P1];
    withinOrBetween=1;
    [outputprime2c.(fields{i}),outputprime2craw.(fields{i})] = HotelingT.HotelingtTest(M, withinOrBetween);
    dob=outputprime2c.(fields{i})(8)*sqrt(2/withinOrBetween);%HotelingTest returns dz for within design
    dtie=0.5*sqrt(2/withinOrBetween);
    k=2;
    n=length(M(:,1));
    lr01prime2c.(fields{i}) = HotelingT.HotelingtlrfromDs(dob,dtie,withinOrBetween,k,n);%HotelingtlrfromDs uses d =t/sqrt(n/2) regardless of the design within or between
    lr01prime2c_cauchy.(fields{i}) = HotelingT.HotelingBF01BasedonD_cauchy(outputprime2c.(fields{i})(8)*sqrt(2/withinOrBetween),k,n,withinOrBetween,1/sqrt(2));
    %lr01= HotelingtlrfromDs(0.25*sqrt(2),0.5*sqrt(2),1,2,46)

    p1=outputprime2c.(fields{i})(7);
    ellipsePara=ellipsePandE20231125(2, 1-p1,P2,'+','k');
    % M=[P2(:,1),P2(:,2)];% one dimension
    % outputprime2cX.(fields{i}) = HotelingtTest(M, 1); %like paried t test
    % [h,p,ci,stats] = ttest2(P2(:,1),P2(:,2));
end
%[outputprime2c,lr01];
%% compare vin with vout and compare vc and vprime to test bias consistency across in and out
T6vin2out= readtable('T6vin2out20240113.csv');
T6vc2prime=readtable('T6vc2prime20240113.csv');
unicateM=table(T6vin2out.Enviornmentname);
[C,ia,ic] = unique(unicateM);
%B=[];
for i=1:length(ia) %((ia) stores unique categoray)
    index=ic==i;%get category i
    key=[C.Var1{i}(1:2)];
    key = erase(key," ");
    %key=key(1:5);
    Vcorrect.(key)=[T6vin2out(index,:).('vin2outcX'), T6vin2out(index,:).('vin2outcY')];
    Vprime.(key)=[T6vin2out(index,:).('vin2outprimeX'), T6vin2out(index,:).('vin2outprimeY')];
    Vinc2prime.(key)=[T6vc2prime(index,:).('vinc2inprimeX'),T6vc2prime(index,:).('vinc2inprimeY')];
    Voutc2prime.(key)=[T6vc2prime(index,:).('voutc2outprimeX'),T6vc2prime(index,:).('voutc2outprimeY')];
    %P.{[]}
    %B.=[B;mean(Oresp(index,:))];       
end

figure
fields=fieldnames(Vprime);
c=jet(length(fields));
% col=zeros(1,3);
p=colormap;
%c=p(int16(256/length(fields))*i+1,:);
%r = randi([1 256],12,1)
for i=1:length(fields)  
    %figure
   
    thiscolor=c(i,:);
    P1=Vcorrect.(fields{i});
    P2=Vprime.(fields{i});
    %scatter(P2(:,1),P2(:,2),['.',c{1}]);
    %scatter(P2(:,1),P2(:,2),'.','ColorVariable',c{1})
    %plot(P2(:,1),P2(:,2),'.','Color',c{1})
    plot(P2(:,1),P2(:,2),'.','Color',thiscolor)
    hold on
    %ellipsePandE20210716(1, 0.95,P2,['+',c{1}])
    output.(fields{i})=ellipsePandE20231125(2, 0.95,P2,'+',thiscolor);
    hold on
    plot(P1(1,1),P1(1,2),'o','Color',thiscolor);
    hold on
    text(P1(1,1),P1(1,2),fields{i},'Color','b');
    Mu=output.(fields{i}).Mu;
    text(Mu(1),Mu(1,2),[fields{i},char(39)],'Color','k');
    %plot(P1(1,1),P1(1,2),['o',c{1}]);
     %do hythothesis test between prime and correct positions using
    %HotelingT
    M=[P2,P1];
    withinOrBetween=1;
    [outputVprime2c.(fields{i}),outputVprime2craw.(fields{i})] = HotelingT.HotelingtTest(M, withinOrBetween);
    
    rawdie=outputprime2craw.([fields{i},'0']).d;
    rawcov=outputVprime2craw.(fields{i}).cov;

    outputVprime2cDie.(fields{i})=sqrt(rawdie*inv(rawcov)*rawdie'); % use the raw distance from the distortion but cov from the consistancy  

    lr01Vprime2c.(fields{i}) = HotelingT.HotelingtlrfromDs(outputVprime2c.(fields{i})(8)*sqrt(2/withinOrBetween),outputprime2c.([fields{i},'0'])(8)*sqrt(2/withinOrBetween),withinOrBetween,k,n);
    lr01Vprime2cUsingRawD.(fields{i}) = HotelingT.HotelingtlrfromDs(outputVprime2c.(fields{i})(8)*sqrt(2/withinOrBetween),outputVprime2cDie.(fields{i})*sqrt(2/withinOrBetween),withinOrBetween,k,n);

    lr01Vprime2c_cauchy.(fields{i}) = HotelingT.HotelingBF01BasedonD_cauchy(outputVprime2c.(fields{i})(8)*sqrt(2/withinOrBetween),k,n,withinOrBetween,1/sqrt(2));
    % p1=outputVprime2c.(fields{i})(3);
    % ellipsePara=ellipsePandE20231125(2, 1-p1,P2,'+','k');
end

figure 
for i=1:length(fields)  
    %figure
   
    thiscolor=c(i,:);
    P1= Vinc2prime.(fields{i});
    P2= Voutc2prime.(fields{i});
    %scatter(P2(:,1),P2(:,2),['.',c{1}]);
    %scatter(P2(:,1),P2(:,2),'.','ColorVariable',c{1})
    %plot(P2(:,1),P2(:,2),'.','Color',c{1})
    plot(P2(:,1),P2(:,2),'.','Color',thiscolor)
    hold on
    %ellipsePandE20210716(1, 0.95,P2,['+',c{1}])
    output.(fields{i})=ellipsePandE20231125(2, 0.95,P2,'+',thiscolor);
    hold on       
    Mu=output.(fields{i}).Mu;
    text(Mu(1),Mu(1,2),[fields{i},'d1'],'Color','k');

    hold on

    output.(fields{i})=ellipsePandE20231125(2, 0.95,P1,'O',thiscolor);
    hold on       
    Mu=output.(fields{i}).Mu;
    text(Mu(1),Mu(1,2),[fields{i},'d0'],'Color','k');

    M=[P2,P1];
    withinOrBetween=1;
    outputVin2out.(fields{i}) = HotelingT.HotelingtTest(M, withinOrBetween);
    lr01Vin2out.(fields{i}) = HotelingT.HotelingtlrfromDs(outputVin2out.(fields{i})(8)*sqrt(2/withinOrBetween),0.5*sqrt(2/withinOrBetween),withinOrBetween,k,n);
    lr01Vin2out_cauchy.(fields{i}) = HotelingT.HotelingBF01BasedonD_cauchy(outputVin2out.(fields{i})(8)*sqrt(2/withinOrBetween),k,n,withinOrBetween,1/sqrt(2));
   
end
 text(0,0,'0,0','Color','r');

 %%
 HotelingT=collectiveHotelingTsquared
 n=(4:1:48)';
 dz=0.5;
 alpha=0.05;
 d=dz*sqrt(2);
 %withinOrBetween=1; 
 k=2;%k is dimension of the data when k is 1 it is similar to power analysis based on d
 withinOrBetween=1;
 %d=0.5;

 [n, df2, critical_F, F, lambdaF, p2]=HotelingT.Hotelingtpowerd(d,k, n,alpha,withinOrBetween);%Hotelingtpowerd uses d =t/sqrt(n/2) regardless of the design within or between

 [n, df2, critical_F, F, lambdaF, p2]

 %% the follow is for test purpose. it shows that two sample hoteling ttest is a special case of manova
 P1= Vinc2prime.('CC')
 P2= Vinc2prime.('Pe');
 M=[P2,P1];
 withinOrBetween=2;
 HotelingT.HotelingtTest(M, withinOrBetween)
 %% save key findings

 [fields12, toutputprime2c]=struct2Matrix(outputprime2c);

 toutput12=array2table(toutputprime2c,'VariableNames',{'Hoteling T','df1','df2','F','Fdf1','Fdf2','p','dz'});

 [toutput12.viewingBuilding, toutput12.B01]=struct2Matrix(lr01prime2c_cauchy);

 toutput12 = movevars(toutput12, 'viewingBuilding', 'Before', 1);


 [fields6, toutputVin2out]=struct2Matrix(outputVin2out);
 toutput6=array2table(toutputVin2out,'VariableNames',{'Hoteling T','df1','df2','F','Fdf1','Fdf2','p','dz'});

 [toutput6.viewingBuilding, toutput6.B01]=struct2Matrix(lr01Vprime2c_cauchy);
 toutput6 = movevars(toutput6, 'viewingBuilding', 'Before', 1);
 

 %save('outputprime2c.mat','toutputprime2c','toutputVin2out','toutput6','toutput12')
 writetable(toutput6,'consistencyTable.xlsx','Sheet','toutput6');
 writetable(toutput12,'consistencyTable.xlsx','Sheet','toutput12');



function [fields, M]=struct2Matrix(outputprime2c)
    M=[]
    fields=fieldnames(outputprime2c);    
    for i=1:length(fields) 
        M= [M;outputprime2c.(fields{i})];
    end
end