
%
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


fields=fieldnames(Pprime);
for i=1:length(fields)  
    
    P1=Pcorrect.(fields{i});
    P2=Pprime.(fields{i});
    
    
    %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)
   
end

 [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);
%[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


fields=fieldnames(Vprime);

for i=1:length(fields)  
   
    P1=Vcorrect.(fields{i}); % correct in to out which should be same for one building across participants
    P2=Vprime.(fields{i}); % estimated in to out
   
     %do hythothesis test between Vprime and Vcorrect 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


[fields6, toutputVprime2c]=struct2Matrix(outputVprime2c);
toutput6=array2table(toutputVprime2c,'VariableNames',{'Hoteling T','df1','df2','F','Fdf1','Fdf2','p','dz'});
[toutput6.viewingBuilding, toutput6.B01]=struct2Matrix(lr01Vprime2c_cauchy);
toutput6 = movevars(toutput6, 'viewingBuilding', 'Before', 1);
toutput6Vprime2c=toutput6;



for i=1:length(fields)  
   
    P1= Vinc2prime.(fields{i}); % in correct to prime
    P2= Voutc2prime.(fields{i}); % out correct to prime
    
 %do hythothesis test between Vinc2prime and using %HotelingT
    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

[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);
toutput6Vin2out=toutput6; 


 %% power analysis
 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)

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