close all;
clear all;
T6vc2prime=readtable('T6vc2prime20240113.csv');
unicateM=table(T6vc2prime.Enviornmentname);
[C,ia,ic] = unique(unicateM);
%B=[];
%%compare prime with correct, testing bias
%tprime2c=table();
keyiCell={};
keyjCell={};
confieldNames={};
for i=1:length(ia) %((ia) stores unique categoray)
    indexi=ic==i;%get category i
    keyi=C.Var1{i}(1:2);
    keyiCell{i}=[keyi,'in'];
    keyjCell{i}=[keyi,'out'];
    for j=1:length(ia)
        indexj=ic==j;
        keyj=C.Var1{j}(1:2);
        keyij=[keyi,keyj];
        xin=T6vc2prime.vinc2inprimeX(indexi);
        yin=T6vc2prime.vinc2inprimeY(indexi);
        xout=T6vc2prime.voutc2outprimeX(indexj);
        yout=T6vc2prime.voutc2outprimeY(indexj);

        dis.(keyij)=sqrt((xin-xout).^2+(yin-yout).^2); %this is the distance or disimilarity of the error vectors
        disAdd.(keyij)=sqrt((xin+xout).^2+(yin+yout).^2); 

        consistency36.(keyij)=(disAdd.(keyij).^2-dis.(keyij).^2)./(disAdd.(keyij).^2+dis.(keyij).^2);
        

        participantID.(keyij)=[T6vc2prime.participant_number(indexi),T6vc2prime.participant_number(indexj)];
        meandis.(keyij)=mean(dis.(keyij)); % this is the mean distance or disimilarity of the error vectors
        meandisAdd.(keyij)=mean(disAdd.(keyij));
        meanconsistency36.(keyij)=mean( consistency36.(keyij));
        consistency36_ttestoutput=myTtestonSample(consistency36.(keyij));
        p1.(keyij)=consistency36_ttestoutput(1);
        t1.(keyij)=consistency36_ttestoutput(2);
        dz1.(keyij)=consistency36_ttestoutput(6);
        % bi_in2out=bi_regressV3Table(xin,yin,xout,yout,...
        % {'d','angle','scale','rsq'});
        % rsqXY.(keyij)= sqrt(bi_in2out.rsq);
        % 
        % rx.(keyij)=corr(xin,xout);
        % ry.(keyij)=corr(yin,yout);

        if strcmp(keyi,keyj)
            confieldNames=[confieldNames,keyij]
        end
    end
end

meandisoutput=inoutMatrixComparison(meandis,keyjCell,keyiCell,'errorSub')
meandisAddoutput=inoutMatrixComparison(meandisAdd,keyjCell,keyiCell,'errorAdd')
meanconsistency36output=inoutMatrixComparison(meanconsistency36,keyjCell,keyiCell,'errorAddsqr-errorSubSqr/errorAddsqr+errorSubSqr')
% dzoutput=inoutMatrixComparison(dz1,keyjCell,keyiCell)
% toutput=inoutMatrixComparison(t1,keyjCell,keyiCell)
% poutput=inoutMatrixComparison(p1,keyjCell,keyiCell)
% rsqXYoutput=inoutMatrixComparison(rsqXY,keyjCell,keyiCell)
% rxoutput=inoutMatrixComparison(rx,keyjCell,keyiCell)
% ry=inoutMatrixComparison(ry,keyjCell,keyiCell)

%confieldNames={'BuBu','CCCC','CaCa','HuHu','PePe','RuRu'} % it could lead to errors if we change names 

[disoutput,disArray]=TestConsistency(dis,confieldNames);
[disAddoutput,disAddArray]=TestConsistency(disAdd,confieldNames);
[consistency36output,consistency36Array]=TestConsistency(consistency36,confieldNames);

consistency36_con_output=myTtestonSample(consistency36Array.dconArray);
consistency36_incon_output=myTtestonSample(consistency36Array.dinconArray);
dconsistency36_con_output=myTtestonSample(consistency36Array.dconArray-consistency36Array.dinconArray);

con_disVsdisAddtestoutput=myTtest(disArray.dconArray,disAddArray.dconArray);
incon_disVsdisAddtestoutput=myTtest(disArray.dinconArray,disAddArray.dinconArray);

consistency_con=(disAddArray.dconArray.^2-disArray.dconArray.^2)./(disAddArray.dconArray.^2+disArray.dconArray.^2);
consistency_incon=(disAddArray.dinconArray.^2-disArray.dinconArray.^2)./(disAddArray.dinconArray.^2+disArray.dinconArray.^2);


consistency_con_ttestoutput=myTtestonSample(consistency_con)
consistency_incon_ttestoutput=myTtestonSample(consistency_incon)
dconsistency_con_ttestoutput=myTtestonSample(consistency_con-consistency_incon)

consistency_conMSE=(disAddArray.conMSE-disArray.conMSE)./(disAddArray.conMSE+disArray.conMSE);
consistency_inconMSE=(disAddArray.inconMSE-disArray.inconMSE)./(disAddArray.inconMSE+disArray.inconMSE);

consistency_conMSE_ttestoutput=myTtestonSample(consistency_conMSE)
consistency_inconMSE_ttestoutput=myTtestonSample(consistency_inconMSE)

dconsistency_conMSE_ttestoutput=myTtestonSample(consistency_conMSE-consistency_inconMSE)

%%
% save consistency score to T2AwithBiregression.csv
T2A=readtable("T2AwithBiregression.csv");
Tcons=struct2table(consistency36Array);
%Tcons=Tcons(:,{'dconArray','dinconArray'}); %change the names to consistencyWithin consistencyBetween
fns=fieldnames(participantID);
TparticipantNo=array2table([participantID.(fns{1})(:,1),consistency36Array.dconArray,consistency36Array.dinconArray],'VariableNames',{'participant_number','consistencyWithin','consistencyBetween'});
TconsParti=[Tcons,TparticipantNo];
T2A=join(T2A, TconsParti);

writetable(T2A,"T2AwithBiregression.csv");
%table([consistency36Array.dconArray,consistency36Array.dinconArray])

%% dconArray=dis.AdAd+dis.CCCC+dis.CaCa+dis.HuHu+dis.PePe+dis.RuRu
function output=myTtest(dconArray,dinconarray)
    [h,p,ci,stats]=ttest(dconArray,dinconarray);
    dz=stats.tstat/sqrt(length(dconArray));
    output=[p,stats.tstat,stats.df,mean(dconArray),mean(dinconarray),dz];
end
function output=myTtestonSample(dconArray)
    [h,p,ci,stats]=ttest(dconArray);
    dz=stats.tstat/sqrt(length(dconArray));
    output=[p,stats.tstat,stats.df,mean(dconArray),0,dz];
end
function [output,arrays]=TestConsistency(dis,confieldNames)
    fields=fieldnames(dis);
    dconArrayTotal=0;
    dinconArrayTotal=0;
    conSSE=0;
    inconSSE=0;
    %confieldNames={'AdAd','CCCC','CaCa','HuHu','PePe','RuRu'};%consistent/within a building
    for i=1:length(fields)
        if ismember(fields(i),confieldNames)
            dconArrayTotal=dconArrayTotal+dis.(fields{i});   
            conSSE=conSSE+dis.(fields{i}).^2
            %dis.AdAd+dis.CCCC+dis.CaCa+dis.HuHu+dis.PePe+dis.RuRu
        else
            dinconArrayTotal=dinconArrayTotal+dis.(fields{i});
            inconSSE=inconSSE+dis.(fields{i}).^2
        end
    end
    
    dconArray=dconArrayTotal/(length(confieldNames));
    dinconarray=dinconArrayTotal/(length(fields)-length(confieldNames));
    [h,p,ci,stats]=ttest(dconArray,dinconarray)
    dz=stats.tstat/sqrt(length(dconArray))
    output=[p,stats.tstat,stats.df,mean(dconArray),mean(dinconarray),dz];
    %[dconArray,dinconarray]
    arrays.dconArray=dconArray;
    arrays.dinconArray=dinconarray;
    arrays.conMSE=conSSE/(length(confieldNames));
    arrays.inconMSE=inconSSE/(length(confieldNames));
    
    % simSqr=1-dconArray.^2./dinconarray.^2;
    % simSqr(simSqr<0)=0;
    % 
    % similarity=sqrt(simSqr);
    % mean(similarity)
    % 
    % sim=1-dconArray./dinconarray;
    % sim(simSqr<0)=0;
    % 
    % %sim=sqrt(simSqr);
    %mean(sim)
end

%% what are the distance of each individual error location to the mean error location of the same building, mean square of which should be the same as mean variance or trace of the error.
% for i=1:length(ia) %((ia) stores unique categoray)
%     indexi=ic==i;%get category i
%     key=C.Var1{i}(1:2);
%     keyin=[key,'in'];
%     keyout=[key,'out'];
% 
%     xin=T6vc2prime.vinc2inprimeX(indexi);
%     yin=T6vc2prime.vinc2inprimeY(indexi);
%     xout=T6vc2prime.voutc2outprimeX(indexi);
%     yout=T6vc2prime.voutc2outprimeY(indexi);
% 
% 
%     % dis2mean.(keyin)=sqrt((xin-mean(xin)).^2+(yin-mean(yin)).^2); %this is the distance or disimilarity to the mean error
%     % dis2mean.(keyout)=sqrt((xout-mean(xout)).^2+(yout-mean(yout)).^2);
% 
%     dis2mean.(keyin)=sqrt((xin-0).^2+(yin-0).^2); %this is the distance or disimilarity to the mean error
%     dis2mean.(keyout)=sqrt((xout-0).^2+(yout-0).^2);
% 
%     % covunbiasedin=mean((xin-mean(xin)).^2+(yin-mean(yin)).^2);
%     % covunbiasedout=mean((xout-mean(xout)).^2+(yout-mean(yout)).^2);
%     covunbiasedin=mean((xin-0).^2+(yin-0).^2);
%     covunbiasedout=mean((xout-0).^2+(yout-0).^2);
% 
%     n=length(xin);
%     varunbiased.(keyin)=covunbiasedin;%cov([xin,yin])*(n-1)/n;
%     varunbiased.(keyout)=covunbiasedout;%cov([xout,yout])*(n-1)/n;   
% 
%     dis.(key)=sqrt((xin-xout).^2+(yin-yout).^2);
% 
%     mse=mean((xin-xout).^2+(yin-yout).^2);
% 
%     rpartialsqr.(keyin)=1-mse/covunbiasedin;
%     rpartialsqr.(keyout)=1-mse/covunbiasedout;
% end
% 
% mean([covunbiased.Adin,covunbiased.Adout,covunbiased.CCin,covunbiased.CCout,covunbiased.Cain,covunbiased.Caout,...
%     covunbiased.Huin,covunbiased.Huout,covunbiased.Pein,covunbiased.Peout,covunbiased.Ruin,covunbiased.Ruout])

% t=struct2table(dis);
% A = table2array(t);
% dinconsistArray=

function [output, h]=inoutMatrixComparison(meandis,keyjCell,keyiCell,title)
    m=stru2matrix(meandis);
    h=distanceMatrix(m,keyjCell,keyiCell,title);
    w=diag(m);
    b= m+diag(diag(m)./diag(m))*1000000;
    b=b(:);
    b(b>10000)=[];
    [h,p,ci,stats]=ttest2(w,b);
    %d=stats.tstat/sqrt(length(w));
    output=[p,stats.tstat,stats.df,mean(w),mean(b)];
end
%distanceMatrix(m,{'Ad','CC','Ca','Hu','Pe','Ru'},{'Ad','CC','Ca','Hu','Pe','Ru'})

function m=stru2matrix(s)
    t=struct2table(s);
    fields=fieldnames(s);    
    total=length(fields);
    rtotal=sqrt(total);
    A = table2array(t);
    m=reshape(A,[rtotal,rtotal]);
end

function T=bi_regressV3Table(X,Y,A,B,VarName)
     meanX=mean(X);
     meanY=mean(Y);
     meanA=mean(A);
     meanB=mean(B);
     res=bi_regressV3(X,Y,A,B,meanX,meanY);
     d=proj39collection.distance([meanX,meanY],[meanA,meanB]);
     
     T=array2table([d,res.regression(5:end)],"VariableNames",VarName);    
end 

function h=distanceMatrix(M,xvalues,yvalues,title)
    figure
    cdata = M;
    %xvalues = {'M1','M2','M3','NoRecovered'};
    %yvalues = {'M1','M2','M3'};
    h = heatmap(xvalues,yvalues,cdata,'ColorLimits',[0 max(cdata,[],"all")]);
    
    % newStr = split(sheetname);
    % h.Title = ['Confusion matrix:',newStr{end-1},' ',newStr{end}];
    % h.XLabel = 'Recovered Model';
    % h.YLabel = 'True Model';
    h.Title =title;
    ax = gca;
    axp = struct(ax);       %you will get a warning
    axp.Axes.XAxisLocation = 'top';
    %savefig([filename(1:end-5),sheetname,'.fig']);
end