% General Purpose:  PCA/(spectral moment approach)

%

% 0. Input:
%   Data should be arranged such that
%
%    RatioData contains all following data/fields
%    name - the name of the study
%    length - the number in the population of the study
%    data  - which contains the data arranged the following way
%                  Investor Data .... Trustee Data
%
%  So RatioData(2).data is the data of the second study: columns label
%                 rounds, rows label which contest
%      each row in data contains [investor play, trustee play]
%   Note: The computer assumes two players,
%       but computes the number of rounds



% 0. (cont)  The user must change the following inputs
%    CreateFigsQ - 1 means to create Figs as the script processes
%                  0 means no Figs
%    IStudyVec   -  an array like [1,3:5] indicating which
%                  studies are to be analyzed
%    RandQ       - use random numbers (1) if yes  ;  Don't change for now
%    SnipSection - such as [2:5]; this allows you to use sections shorter
%                    than the three rounds.
%                  Default is [1:6]; Don't change for now
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%  
% Start in cd  W:\prb\Brooks\RandomControls\Oct062007
%


cd W:\prb\Brooks\Shuffled

load Sept10th2007.mat RatioData NameForStudyVec



CreateFigsQ=0; % Make 0 if you don't want figs; 1 if you do want
IStudyVec=[1,3:16];
% IStudyVec=[9:10];
RandQ=0;           % Don't change for now
SnipSection=[1:6]; % Don't change for now
NumLoops = 100;  %


%   Variables created:
%
%   NameForStudyVec    the list of names describing each study
%   NPeopleVec         (the number of contests in each population)
%   NumPeople          (the total number of contests)
%   NumRounds          may not necessarily be 10
%   RatioDataNMDS      (all of the investor/trustee data)
%   Snips               a matrix of all the three round
%   SnipsperTrustee     averaged to level of Trustee
%   SnipsperContest     a matrix of all the three round
%   Fisher12            is the value of the Fisher Discriminant taken
%                       pairwise between populations
%
%   There is a lot more things, but a lot of it is self explanatory


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%         Begin Section 1  Calculate Snips, SnipsPerTrustee, SnipsPerStudy
%           NPeopleVec, ImeanSnips, TmeanSnips, IstdSnips, TstdSnips


NumLoops= (RandQ>0)*NumLoops; % This is really the number of extra loops

NumStudies =length(IStudyVec);
NameForStudyVec = cell(1,NumStudies);

clear Fisher12Monster
for jLoop=1:(NumLoops+1)
    jLoop
    countPeople=0;
    for jStudy=1:NumStudies
        jStudyIndex=IStudyVec(jStudy);
        NPeople =RatioData(jStudyIndex).length;
        NPeopleVec(jStudy)=NPeople;
        temp = RatioData(jStudyIndex).data;
        tempSize= size(temp);
        if ((RandQ==1)& (NumLoops+1-jLoop)), temp=rand(tempSize); end
        RatioDataNMDS( (countPeople+1):(countPeople+NPeople),1:20) = temp;
        countPeople=countPeople+NPeople;
        NameForStudyVec{jStudy} = RatioData(jStudyIndex).name;
    end

    NumPeople =countPeople;

    sizeRatioData = size(RatioDataNMDS);
    NumRounds = sizeRatioData(2)/2;
    NumSnipsPerContest=NumRounds-2;

    clear countPeople NPeople jStudy jStudyIndex temp tempSize sizeRatioData

    Snips= zeros(NumPeople*8,6);

    for jPeople=1:NumPeople
        for jRound=3:10
            Snips( (jPeople-1)*8+jRound-2, 1:2:5) = RatioDataNMDS(jPeople,(jRound-2):jRound );
            Snips( (jPeople-1)*8+jRound-2, 2:2:6) = RatioDataNMDS(jPeople,(jRound+8):(jRound+10));
        end
    end


    NumSnips = NumPeople*8;

    temp=Snips;
    clear Snips;
    for jSnips=1:NumSnips
        Snips(jSnips,:) = temp(jSnips,SnipSection);
    end


    clear SnipsPerStudy SnipsPerTrustee
    for jStudy=1:NumStudies
        NPeople=NPeopleVec(jStudy);
        FirstPerson=  sum(NPeopleVec(1:(jStudy-1)))+1;
        LastPerson = FirstPerson-1+NPeople;
        for jPeople=FirstPerson:LastPerson
            SnipsPerTrustee(jPeople,:) = mean( Snips( ((jPeople-1)*NumSnipsPerContest+1): ...
                (jPeople*NumSnipsPerContest),:));
        end
        SnipsPerStudy(jStudy,:) = mean( SnipsPerTrustee( FirstPerson:LastPerson ,:));
    end




    clear jSnips jPeople jRound  temp jStudy FirstPerson LastPerson

    LSnipSection = length(SnipSection);

    if (LSnipSection==2), IIVec=[1]; TTVec=[2]; end
    if (LSnipSection==4), IIVec=[1,3]; TTVec=[2,4]; end
    if (LSnipSection==6), IIVec=[1,3,5]; TTVec=[2,4,6]; end


    facTor =  length(IIVec)*sqrt(NumSnips);

    ImeanSnips = mean(mean(Snips(:,IIVec)));
    IstdSnips =0;
    for jj=1:length(IIVec)
        IstdSnips = IstdSnips +  norm(Snips(:,IIVec(jj ) ) - ImeanSnips)/facTor ;
    end

    TmeanSnips = mean(mean(Snips(:,TTVec)));
    TstdSnips = 0;
    for jj=1:length(IIVec)
        TstdSnips = TstdSnips +  norm(Snips(:,TTVec(jj ) ) - TmeanSnips)/facTor ;
    end


    clear facTor jj

clear zscoreSnips  zscoreSnipsPerStudy zscoreSnipsPerTrustee

    zscoreSnips(:,IIVec) = (Snips(:,IIVec) - ImeanSnips)/ IstdSnips  ;
    zscoreSnips(:,TTVec) = (Snips(:,TTVec) - TmeanSnips)/ TstdSnips  ;

    zscoreSnipsPerStudy(:,IIVec) = (SnipsPerStudy(:,IIVec) - ImeanSnips)/ IstdSnips  ;
    zscoreSnipsPerStudy(:,TTVec) = (SnipsPerStudy(:,TTVec) - TmeanSnips)/ TstdSnips  ;

    zscoreSnipsPerTrustee(:,IIVec) = (SnipsPerTrustee(:,IIVec) - ImeanSnips)/ IstdSnips  ;
    zscoreSnipsPerTrustee(:,TTVec) = (SnipsPerTrustee(:,TTVec) - TmeanSnips)/ TstdSnips  ;


    zscoreSnips = Snips;





    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% end Section 1%%%%%%%%%%%%%%%%%%%%%

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %                         Begin Section 2. PCA
    %

    % -----------------------------Begin Section 2 PCA

    [coefsSnips,latentSnips,explainedSnips] =pcacov(zscoreSnips'*zscoreSnips);
    latentSnips= round(latentSnips); % []
    %  explainedSnips
    signsSnips= sign(mean(zscoreSnips*coefsSnips));
    coefsSnips = coefsSnips * diag(signsSnips);

    if CreateFigsQ
        saveasString= 'SnipEigenVectors';
        figure('name',saveasString);
        colormap(gray);
        bar(coefsSnips');
        set(gca,'FontName','Arial','FontSize',18,'FontWeight','bold','box','on')
        title(saveasString);
        explainedtemp = round(10*explainedSnips)/10;
        set(gca,'XTickLabel', {num2str(explainedtemp)})
        set(gca,'YTick', [-1,-0.8,-.6,-0.4,-.2,0.0,.2,0.4,.6,0.8,1.])
        set(gca,'LineWidth',3)
        set(gca,'YTickLabel', {'','-0.8','','-0.4','','0.0','','0.4','','0.8',''})
        saveas(gcf,strcat(saveasString,'.jpg'),'jpg');
        saveas(gcf,strcat(saveasString,'.fig'),'fig') ; saveas(gcf,strcat(saveasString,'.eps'),'eps');
    end

    [coefsGames,latentGames,explainedGames] =pcacov(zscoreSnipsPerTrustee'*zscoreSnipsPerTrustee);
    explainedGames = round(100*explainedGames)/100;
    %  [85.57;12.07;1.6;0.48;0.16;0.11;]
    signsGames= sign(mean(zscoreSnipsPerStudy*coefsGames));
    coefsGames = coefsGames * diag(signsGames);

    if CreateFigsQ
        saveasString= 'Game EigenVectors';
        figure('name',saveasString);
        colormap(gray); bar(coefsGames');
        set(gca,'FontName','Arial','FontSize',22,'FontWeight','bold','box','on')
        title(saveasString);
        explainedtemp = round(explainedGames); set(gca,'XTickLabel', {num2str(explainedtemp)})
        set(gca,'YTick', [-1,-0.8,-.6,-0.4,-.2,0.0,.2,0.4,.6,0.8,1.])
        set(gca,'LineWidth',3,'YTickLabel', {'','-0.8','','-0.4','','0.0','','0.4','','0.8',''})
        saveas(gcf,strcat(saveasString,'.fig'),'fig') ; saveas(gcf,strcat(saveasString,'.eps'),'eps');
    end

    scoresS = zscoreSnips * coefsSnips;
    scoresS2 = scoresS'*scoresS;

    scoresG = zscoreSnipsPerTrustee * coefsSnips;
    scoresG2 = scoresG'*scoresG;


    %-----------------------------------------

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %
    %----------------  Section 2 continued:  Construct symmetric eigenvectors
    e1  = -ones(1,LSnipSection);
    e1 = e1/norm(e1);
    %
    scores1S = zscoreSnips* e1';
    scores1S'*scores1S/trace(scoresS2)  % gives 56.59,  16.58% for random vectors 07/25/07
    latent1 = scores1S'*scores1S;%
    scores1G = zscoreSnipsPerTrustee*e1';
    scores1G'*scores1G/trace(scoresG2)  % gives 84.95, 45.89% for random vectors

    %-----------------------------------------
    e2  = e1;
    for jj=2:2:LSnipSection , e2(jj)=- e1(jj); end
    %
    scores2S = zscoreSnips* e2';
    scores2S'*scores2S/trace(scoresS2)  % gives 11.36,  16.66% for random vectors
    latent2 = scores2S'*scores2S;
    scores2G = zscoreSnipsPerTrustee*e2';
    scores2G'*scores2G/trace(scoresG2)  % gives 12.53, 43.07% for random vectors



    zscoreSnipsWO2= zscoreSnips- zscoreSnips*e1'*e1 -zscoreSnips*e2'*e2;


    [coefsSnipsWO2,latentSnipsWO2,explainedSnipsWO2] =pcacov(zscoreSnipsWO2'*zscoreSnipsWO2);

    % latentSnipsWO2 2040.4, 1809.3, 1259.3, 1157.2  for shuffled vectors
    % latentSnipsWO2 2124.4, 2072.3, 1709.3, 1627.2  for shuffled all vectors
    % latentSnipsWO2 3406     3314     3200     3117 for random vectors
    latentNormSnips = latent1+latent2+ sum(latentSnipsWO2 );
    % latentNormSnips 9585.4, 19,535 for random vectors

    clear temp jStudy jStudyIndex countPeople
%   ----------------------------------------------------------------
%  Do Canonical correlation analysis at this point;   
%   Large eigenvector is investor reciprocity
   
    InvSnipsWO2= zscoreSnipsWO2(:,1:2:6);
    TruSnipsWO2= zscoreSnipsWO2(:,2:2:6);

    
    IMat= InvSnipsWO2'*InvSnipsWO2;
    TMat= TruSnipsWO2'*TruSnipsWO2;
    ZMat= InvSnipsWO2'*TruSnipsWO2;
    
    et1= [-1 0 1 ] /sqrt(2);
    et2= [1 -2 1 ] /sqrt(6);
    
    IMatp=zeros(2,2); TMatp=zeros(2,2); ZMatp=zeros(2,2);
    IMatp(1,1)=et1*IMat*et1';  IMatp(1,2)=et1*IMat*et2';
    IMatp(2,1)=et2*IMat*et1';  IMatp(2,2)=et2*IMat*et2';
    TMatp(1,1)=et1*TMat*et1';  TMatp(1,2)=et1*TMat*et2';
    TMatp(2,1)=et2*TMat*et1';  TMatp(2,2)=et2*TMat*et2';
    ZMatp(1,1)=et1*ZMat*et1';  ZMatp(1,2)=et1*ZMat*et2';
    ZMatp(2,1)=et2*ZMat*et1';  ZMatp(2,2)=et2*ZMat*et2';

    
    IMatH = IMatp^(-.5);
    TMatH = TMatp^(-.5);

    ZHat = IMatH*ZMatp*TMatH;
  
    ZHatI = ZHat *ZHat';     [WI,DI] = eig(ZHatI)
    ZHatT = ZHat'*ZHat ;     [WT,DT] = eig(ZHatT)
    
    vIpmax = IMatH*WI(:,2);
    vTpmax = TMatH*WT(:,2);
    
    vImax= et1*vIpmax(1) +et2*vIpmax(2); vImax=vImax/norm(vImax)
    vTmax= et1*vTpmax(1) +et2*vTpmax(2); vTmax=vTmax/norm(vTmax)
    vIT(1:2:6)=vImax;  vIT(2:2:6) = vTmax;
    figure
    
    barRange=1:2:5; bar(1:2:5,vIT(barRange))
    hold on;
    barRange=2:2:6; bar(barRange,vIT(barRange))
%    bar(vIT)
    AxelIR =  (vImax(2)-vImax(3))*(vTmax(1)-vTmax(2))
    AxelTR = ((vImax(1)-vImax(2))*(vTmax(1)-vTmax(2)) +  ...
              (vImax(1)-vImax(2))*(vTmax(1)-vTmax(2))  )/2

    %-----------------------------------------
    % first 4 eigenvectors describe 87%
    %

    e3 = coefsSnipsWO2(:,1); e3=e3';
    e4 = coefsSnipsWO2(:,2); e4=e4';
    e5 = coefsSnipsWO2(:,3); e5=e5';
    e6 = coefsSnipsWO2(:,4); e6=e6';


    %%%%           From here out use full 3 rounds

    Projection= [e1',e2',e3',e4',e5',e6'];


    signsGames= sign(mean(zscoreSnips * Projection));
    Projection = Projection * diag(signsGames);


    zscoreSnipsTot = zscoreSnips * Projection;

    zscoreSnipsTot2 = zscoreSnipsTot'  * zscoreSnipsTot ;

    tzscoreSnipsTot2 = trace(zscoreSnipsTot2);
    for j=1:6
        expSnipsTot2(j) =  zscoreSnipsTot2(j,j)/ tzscoreSnipsTot2;
    end


    clear j

    if CreateFigsQ
        saveasString= 'FinalEigenvectors' ;%open(strcat(saveasString,'.fig')) ;
        figure('name',saveasString);
        colormap(gray);
        bar(Projection');
        set(gca,'FontName','Arial','FontSize',18,'FontWeight','bold','box','on')
        set(gca,'LineWidth',3)
        title(saveasString);
        explainedtemp = round(1000*expSnipsTot2)/10;
        set(gca,'YLim',[-1 1]);
        set(gca,'XTickLabel', {num2str(explainedtemp' )})
        set(gca,'YTick', [-1,-0.8,-.6,-0.4,-.2,0.0,.2,0.4,.6,0.8,1.])
        set(gca,'YTickLabel', {'','-0.8','','-0.4','','0.0','','0.4','','0.8',''})
        saveas(gcf,strcat(saveasString,'.fig'),'fig') ;
        saveas(gcf,strcat(saveasString,'.eps'),'eps');
        saveas(gcf,strcat(saveasString,'.jpg'),'jpg');
    end
    %


    Projection(:,3)= vIT';


    %-----------------------------------------

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %
    %         Section 3: Build Classification Vector

    % reassemble zscoreSnipsTot into InvTruData

    for j=1:NumPeople
        for jRound=1:8
            rHi=jRound*6;
            rLo=rHi-5;;
            InvTruData(j,(rLo:rHi))=zscoreSnips(8*(j-1)+jRound,1:6);
        end
    end, clear j jRound rHi rLo


    for j=1:NumPeople
        for jmoment =1:6
            x=InvTruData(j,:);
            moment=jmoment-1;
            InvTruCum(j,jmoment) =  CumDiscrete(x,Projection,jmoment);
        end
        RecipInvVec(j)    = Reciprocity(x,6,0) ;
        RecipTruVec(j)    = Reciprocity(x,6,1) ;
        ResponseTruVec(j) = Reciprocity(x,6,2) ;
        ResponseInvVec(j) = Reciprocity(x,6,3) ;
        RegretVec(j)      = Reciprocity(x,6,6) ;
        EnvyVec(j)        = Reciprocity(x,6,7) ;
        GuiltVec(j)       = Reciprocity(x,6,8) ;
    end, clear x j jmoment moment

    RecipAveVec= (RecipInvVec+ RecipTruVec)/2;
    RecipDffVec= (RecipInvVec -RecipTruVec)/2;

     clear ClassVec
     ClassVec=InvTruCum(:,1:3);     % good
     ClassVec(:,4)=RegretVec(:);  % good
     ClassVec(:,5)=EnvyVec(:);  % good
%      ClassVec(:,1:2)=[];  % Try not using means
%      ClassVec(:,2:3)=[];  % Try not using Regret Envy
    
    temp=size(ClassVec); 
    LClassVec=temp(2); clear temp
    


    clear j  jj jRound rHi rLo jmoment


    mClassVec= mean(ClassVec);
    sClassVec=std(ClassVec);

    clear zClassVec
    for j=1:NumPeople
        zClassVec(j,:) =  (ClassVec(j,:)  - mClassVec)./sClassVec;
    end , clear j


    clear  Contest sContest semContest signifContest % average over people to get contest data
    for jStudy=1:NumStudies
        NPeople     = NPeopleVec(jStudy);
        FirstPerson = sum(NPeopleVec(1:(jStudy-1)))+1;
        LastPerson  = FirstPerson+ NPeople -1;
        Contest(jStudy,:)  = mean(zClassVec(  FirstPerson:LastPerson,:));
        sContest(jStudy,:) = std(zClassVec(   FirstPerson:LastPerson,:));
        semContest(jStudy,:) = sContest(jStudy,:) /sqrt(NPeople);
        signifContest(jStudy,:) =   Contest(jStudy,:)./  semContest(jStudy,:);
        Contest(jStudy,:)  = mean(zClassVec(  FirstPerson:LastPerson,:));
   end;   clear jStudy NPeople FirstPerson LastPerson

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %         Section 4: Do Linear Discriminant Analysis



    clear Fisher12 PerCentSplit AxesMat

    
    jjPairs = [2 13;13 15; 8 13;9 13; 4 13;13 14; 8 9];
    jjPairs = [8 13 ];
    sjjPairs=size(jjPairs);
    jjjkVec=[];
    for jjP=1:sjjPairs(1)
        jj=jjPairs(jjP,1);
        jk=jjPairs(jjP,2);
        jjjkVec = [ jjjkVec ((jj-1)*NumStudies+jk)];
    end
        
    axestot = zeros(1,LClassVec);  % change me
    for jj=1:NumStudies
        NP1  = NPeopleVec(jj);
        FP1  = sum(NPeopleVec(1:(jj-1)))+1;
        LP1  = FP1+ NP1 -1;
        Fisher12(jj,jj) = 0;
        for jk=(jj+1):NumStudies
            NP2  = NPeopleVec(jk);
            FP2  = sum(NPeopleVec(1:(jk-1)))+1;
            LP2  = FP2+ NP2 -1;
            NPVec= [NP1 NP2] ;
            clear Block
            Block(1:NP1,:)              = zClassVec(FP1:LP1,:)  ;
            Block((NP1+1):(NP1+NP2), :) = zClassVec(FP2:LP2,:) ;
            [ axes,signif, separations ] = LDA( Block, NPVec);
%            axes=[0 0 1 0]';  % CHANGEME 9/27/2007
            BlockProj  = Block* axes;
            BlockProj1 = BlockProj(1:NP1);
            BlockProj2 = BlockProj((NP1+1):(NP1+NP2));
            [n1,xout1]=  hist(BlockProj1);
            [n2,xout2]=  hist(BlockProj2);
            n1Tot=sum(n1);
            n2Tot=sum(n2);
            xoutTot=sort([xout1,xout2]);
            nxoutTot=length(xoutTot);
            NumSplit12Max=0;
            for jz=1:(nxoutTot-1)
                split=(xoutTot(jz)+xoutTot(jz+1))/2;
                n1Less= length( find(BlockProj1<split));
                n1More= length( find(BlockProj1>split));
                n2Less= length( find(BlockProj2<split));
                n2More= length( find(BlockProj2>split));
                vv = [ n1Less/n1Tot+n2More/n2Tot,n2Less/n2Tot+n1More/n1Tot];
                [NumSplit12,Indices] = max(vv);
                if (NumSplit12>NumSplit12Max)
                    NumSplit12Max=NumSplit12;
                    vvmax=vv;
                    IndexMax=Indices(1);
                    nsMax= [ n1Less n1More n2Less n2More];
                end
            end
            PerCentSplit(jj,jk)=NumSplit12Max*50;% /sum(nsMax);
            PerCentSplit(jk,jj)= PerCentSplit(jj,jk);
            % Break up BlockProj into 100 pieces. Find where it
            if (sum(((jj-1)*NumStudies+jk)==jjjkVec)  & jLoop>NumLoops) % (2,13), (5,9) (8,13)   %   CHANGE ME
                if 1 % CreateFigsQ
                    titleString = strcat(NameForStudyVec(jj),': vs :',NameForStudyVec(jk));
                    saveString  = strcat(num2str(jj),'_',num2str(jk) );
                    figure('name',saveString);
                    n1=100*n1/NP1
                    n2=100*n2/NP2
                    bar(xout1,n1,'BarWidth',0.2); hold
                    bar(xout2,n2,'FaceColor','r','BarWidth',0.2); hold
                    title(titleString,'FontSize',14,'FontWeight','bold');
                    set(gca,'FontName','Arial','FontSize',16,'FontWeight','bold','box','on')
                    set(gca,'LineWidth',3)
                    %            set(gca,'YLim',[-1 1]);
                    %            set(gca,'XTickLabel', {num2str(explainedtemp)})
                    %            set(gca,'YTick', [-1,-0.8,-.6,-0.4,-.2,0.0,.2,0.4,.6,0.8,1.])
                    %            set(gca,'YTickLabel', {'','-0.8','','-0.4','','0.0','','0.4','','0.8',''})
                    saveas(gcf,strcat(saveString,'.fig'),'fig') ;
                    saveas(gcf,strcat(saveString,'.jpg'),'jpg')
                    % set(get(h(1),'BaseLine'),'LineWidth',2,'LineStyle',':')
                    axes
                   pause
                end
            end
            Fisher12(jj,jk) = signif;
            Fisher12(jk,jj) = signif;
            axestot= axestot+abs(axes');
%
            Inds(jj,jk)= .5*(2*NumStudies+1-jj)*jj-NumStudies+(jk-jj);
            InvIndsjj(   .5*(2*NumStudies+1-jj)*jj-NumStudies+(jk-jj))= jj;
            InvIndsjk(   .5*(2*NumStudies+1-jj)*jj-NumStudies+(jk-jj))= jk;
            AxesMat(Inds(jj,jk),:) = axes';
            if (NumLoops+1-jLoop)
                Fisher12Monster(jLoop,Inds(jj,jk)) = signif;
            end % 
        end     % jk
    end         % jj

end %   end Monster Loop (for noise trials)
clear jj jStudy NPeople jPeople jZ FP1 FP2 LP1 LP2 n1 n2 jk jz



%% ______________________
%   Wrote this for Read   Feb 2008


range= -3.2:0.4:3.2; % for Per Imp
range= -2.4:0.3:2.4; % for BPD Imp
hold off;
hist(BlockProj1,range);
hold on;
hist(BlockProj2,range);
cc1=[.655 .867 1.];  % BPD           These are colors for the different populations
cc2=[.799 .597 1.];   % Imp
cc3=[.358 .596 1.]; cc2=cc3

childs= get(gca,'Children');
FirstChild = childs(1);
outputX=get(FirstChild,'XData');
outputY=get(FirstChild,'YData');
outputY = outputY*200/sum(sum(outputY));
set(FirstChild,'FaceAlpha',0.5);      %  These are transparencies
set(FirstChild,'FaceColor',cc1/norm(cc1));
set(FirstChild,'XData',outputX)
set(FirstChild,'YData',outputY)

SecondChild = childs(2);  % This is Imp both times
outputX=get(SecondChild,'XData');
outputY=get(SecondChild,'YData');
outputY = outputY*200/sum(sum(outputY));
set(SecondChild,'FaceAlpha',0.5);
set(SecondChild,'FaceColor',cc2/norm(cc2));
set(SecondChild,'XData',outputX); % after changing anything, one has to set the correct outputX etc
set(SecondChild,'YData',outputY); %    back to the plot
title(titleString,'FontSize',14,'FontWeight','bold');
set(gca,'FontName','Arial','FontSize',16,'FontWeight','bold','box','on')
set(gca,'LineWidth',3)

legend('BPD-M', 'Imp')
legend('Per', 'Imp')

legend('boxoff')




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%




% Fisher12=Fisher12MainB;


ProbJustRandom=zeros(NumStudies);

for jj=1:NumStudies  %  calculate cdf
    for jk=(jj+1):NumStudies
        tempRand= Fisher12Monster(:,Inds(jj,jk)) ;
        stempRand=sort(tempRand);
        val = Fisher12(jj,jk);
        NumLarger= length(find(stempRand>val));
        ProbJustRandom(jj,jk)=100*NumLarger/NumLoops;
        ProbJustRandom(jk,jj)=100*NumLarger/NumLoops;
    end
end

j=9; range1= [ (1+sum(NPeopleVec(1:(j-1)))):(sum(NPeopleVec(1:j)))]
j=13;range2= [ (1+sum(NPeopleVec(1:(j-1)))):(sum(NPeopleVec(1:j)))]
zz1= zClassVec(range1,:);
zz2= zClassVec(range2,:);

figure; 
subplot(1,2,1); imagesc(zz1); colormap(gray);
subplot(1,2,2); imagesc(zz2); colormap(gray);








%             Bootstrapping experiment







clear Fisher12p
axestot = zeros(1,LClassVec);  % change me

for jj=1:NumStudies
    jj
    %            jj=8;
    NP1  = NPeopleVec(jj);
    FP1  = sum(NPeopleVec(1:(jj-1)))+1;
    LP1  = FP1+ NP1 -1;
    for jk=(jj+1):NumStudies
        %                jk=13;
        NP2  = NPeopleVec(jk);
        FP2  = sum(NPeopleVec(1:(jk-1)))+1;
        LP2  = FP2+ NP2 -1;
        NPVec= [NP1 NP2] ;
        clear Block
        BlockI(1:NP1,:)              = ClassVec(FP1:LP1,:)  ;
        BlockI((NP1+1):(NP1+NP2), :) = ClassVec(FP2:LP2,:) ;
        %           shuffle the last 3 columns of Block

        for jBoot=1:1000
            vv3 = randperm(NP1+NP2);
            vv4 = randperm(NP1+NP2);
            vv5 = randperm(NP1+NP2);
            %
            Block3= BlockI(vv3,3);
            Block4= BlockI(vv3,4);
            Block5= BlockI(vv3,5);
            %
            Block(:,3) =Block3;
            Block(:,4) =Block4;
            Block(:,5) =Block5;
            %
            [ axes,signif, separations ] = LDA( Block, NPVec);
            %            axes=[0 0 1 0]';  % CHANGEME 9/27/2007
            BlockProj  = Block           * axes;
            BlockProj1 = BlockProj(1:NP1)      ;
            BlockProj2 = BlockProj((NP1+1):(NP1+NP2));
            [n1,xout1]=  hist(BlockProj1)      ;
            [n2,xout2]=  hist(BlockProj2)      ;
            n1Tot=sum(n1);
            n2Tot=sum(n2);
            xoutTot=sort([xout1,xout2]);
            nxoutTot=length(xoutTot);
            if (0& jj==2 & jk==13 & jLoop>NumLoops) % (2,13), (5,9) (8,13)   %   CHANGE ME
                if 1 % CreateFigsQ
                    titleString = strcat(NameForStudyVec(jj),': vs :',NameForStudyVec(jk));
                    saveString  = strcat(num2str(jj),'_',num2str(jk) );
                    figure('name',saveString);
                    n1=100*n1/NP1
                    n2=100*n2/NP2
                    bar(xout1,n1,'BarWidth',0.2); hold
                    bar(xout2,n2,'FaceColor','r','BarWidth',0.2); hold
                    title(titleString,'FontSize',14,'FontWeight','bold');
                    set(gca,'FontName','Arial','FontSize',16,'FontWeight','bold','box','on')
                    set(gca,'LineWidth',3)
                    %            set(gca,'YLim',[-1 1]);
                    %            set(gca,'XTickLabel', {num2str(explainedtemp)})
                    %            set(gca,'YTick', [-1,-0.8,-.6,-0.4,-.2,0.0,.2,0.4,.6,0.8,1.])
                    %            set(gca,'YTickLabel', {'','-0.8','','-0.4','','0.0','','0.4','','0.8',''})
                    saveas(gcf,strcat(saveString,'.fig'),'fig') ;
                    saveas(gcf,strcat(saveString,'.jpg'),'jpg')
                    % set(get(h(1),'BaseLine'),'LineWidth',2,'LineStyle',':')
                    axes
                    %                   pause
                end
            end
            %            Fisher12(jj,jk) = signif;
            %            Fisher12(jk,jj) = signif;
            axestot= axestot+abs(axes');
            %
            Inds(jj,jk)= .5*(2*NumStudies+1-jj)*jj-NumStudies+(jk-jj);
            InvIndsjj(   .5*(2*NumStudies+1-jj)*jj-NumStudies+(jk-jj))= jj;
            InvIndsjk(   .5*(2*NumStudies+1-jj)*jj-NumStudies+(jk-jj))= jk;
            AxesMat(Inds(jj,jk),:) = axes';
            Fisher12Boot(jBoot) = signif;
            Fisher12p(jj,jk) =length(find(Fisher12Boot>Fisher12(jj,jk))) ;
            Fisher12p(jk,jj) = Fisher12p(jj,jk);

        end % jBoot
    end     % jk
end         % jj



Fisher12Main= Fisher12;

% Fisher12=Fisher12MainB;

Fisher12Red=round(Fisher12*100 )/100
NameForStudyVecRed = NameForStudyVec;
PerCentSplitRed =  PerCentSplit;
Fisher12pRed=Fisher12p; 

for jj=NumStudies:(-1):1
    if (NPeopleVec(jj)<10)
        Fisher12Red(:,jj) = [];
        Fisher12Red(jj,:) = [];
        Fisher12pRed(:,jj) = [];
        Fisher12pRed(jj,:) = [];
        NameForStudyVecRed(jj)=[];
        PerCentSplitRed(:,jj)=[];
        PerCentSplitRed(jj,:) = [];
    end
end


Fisher12Red = round(20*Fisher12Red)/10;
lFisher12Red =  length(Fisher12Red);

vjk  = find(Fisher12Red>1)   ; lvjk=length(vjk);
wjk  = find(Fisher12pRed<50); lwjk=length(wjk);
vTot = find(Fisher12Red<0); 


Count=0; vwAccept=[]; vwReject=[];
for jj=1:lvjk
    vj =   1+ mod(vjk(jj)-1, lFisher12Red);
    vk= 1+ (vjk(jj)-vj)/lFisher12Red;
    if vj<vk
        vwjkInd= wjk(find(wjk==vjk(jj)));
        vwjkInd
        if length(vwjkInd)
            Count=Count+1;
            [Count,vj,vk];
            NameForStudyVecRed([vj vk]);
            vwAccept= [vwAccept vjk(jj)];
        end
    end
end

% I find this splits 22 (Count) pairs out of a possible of 
%    10*9/2 =45
    
        

% What numbers appear in same list



 
[ positions, Delta ] = nmds_ui( Fisher12Red, 100 );

saveasString= 'NMDS' ;
figure('name',saveasString);
scatter(positions(1,:),positions(2,:));
gname(NameForStudyVecRed)
colormap(gray);
set(gca,'FontName','Arial','FontSize',22,'FontWeight','bold','box','on')
set(gca,'LineWidth',3)
title(saveasString);
saveas(gcf,strcat(saveasString,'.fig'),'fig') ;
saveas(gcf,strcat(saveasString,'.eps'),'eps')

 
Kendall =  computeKendall(positions, Fisher12Red);
 
Ncorrect= sqrt(2* sqrt(2 *max([Kendall,0])))

% Ncorrect =    10.8776
% Kendall  =    1750


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%