%patterns: probabilities and expectation
%choose d=number of bits, p=probability of 1
%first table: expected waiting times
%second table: probability P(iBj) that pattern i 
%  occurs before pattern j
%third table: all tripples i,j,k such that 
%  P(iBj), P(iBk)<1/2, but i wins in the three-way game, 
%  that is:
%  P(iBjVk)>P(jBiVk)>P(kBiVj),
%  where P(iBjVk) is the prob. that the pattern 
%  i occurs before either pattern j or pattern k
%  label *: also, P(jBk)<0.5 but P(jBiVk)>P(kBiVj), 
%  (or vice versa) so that the order in the three-way 
%  game is reversed compared to the two-way game.

global p q;
p=1/2;
d=2; 


q=1-p;
 
% z=exp1([0 1 0 0],[1 0 0 1]);
% z=pf([1 0 0 1],[0 1 0 0]);
% z=pf([1 0 1 0],[0 1 0 0]);
% z=exp0([0 0 0]);
%z=exp0([1 0 1 1 1 0 1])
%z=pf3( [1 0 1 1],[0 1 0 1], [1 0 0 0]);
%z=pf([0 0 0 0],[1 1 1 1 0])



fprintf('---waiting expectations:\n'); 
for i=0:(2^d-1)
   a=dec2bin(i);
   a=str2num(a(:))'; 
   na=size(a); na=na(2); 
   a=[zeros(1,(d-na)) a]; 
   fprintf('%2d:',i); 
   fprintf('%1.4f\n',exp0(a));  
end

fprintf('---one-on-one probabilities:\n');

cl=zeros(2^d);

 for j= 0:(2^d-1)
     fprintf('     %2d', j); 
 end
  fprintf('\n'); 
for i=0:(2^d-1)
    fprintf('%2d:',i); 
    for j= 0:(2^d-1)
        a=dec2bin(i);
        a=str2num(a(:))'; 
        na=size(a); na=na(2); 
        a=[zeros(1,(d-na)) a]; 
        b=dec2bin(j);
        b=str2num(b(:))'; 
        nb=size(b); nb=nb(2);
        b=[zeros(1,(d-nb)) b]; 
%         s=dec2bin(bin2dec(num2str(a))+bin2dec(num2str(b)))
%         out=str2num(s(:))
%        fprintf('%1.4f ',pf(a,b)); 
        eps=0.0000001;
        if ((pf(a,b)>0.5+eps)&& (exp0(a)>exp0(b)))
            cl(i+1,j+1)=1;
        end
        if ((pf(a,b)<0.5-eps)&& (exp0(a)<exp0(b)))
            cl(i+1,j+1)=1;
        end
        if (cl(i+1,j+1)==1) 
            fprintf('%1.4f*',pf(a,b)); 
        else
            fprintf('%1.4f ',pf(a,b)); 
        end
    end
    fprintf('\n');
end

z=pf([1 0 0],[0 0 1]);
z=pf3([0 1 0 1], [1 1 0 1], [1 1 0 0]);
z=pf3([1 0 0 0], [0 1 0 1],  [1 0 1 1]);
z=pf3( [1 0 1 1],[0 1 0 1], [1 0 0 0]);

% notrials=10000; 
% for tr=1:notrials
%     i=randi([0,2^d-1]);
%     j=randi([0,2^d-1]);
%     k=randi([0,2^d-1]);

fprintf('---voting paradox triples:\n');
fprintf(' i  j  k  P(iBj)   P(iBk)   P(jBk)  | P(iBjVk) P(jBiVk) P(kBiVj)\n');

for (i=1:2^d-1)
for (j=1:2^d-1)
for (k=1:2^d-1)
if ((i~=j)&&(i~=k)&&(j<k))
    a=dec2bin(i);
        a=str2num(a(:))'; 
        na=size(a); na=na(2); 
        a=[zeros(1,(d-na)) a]; 
    b=dec2bin(j);
        b=str2num(b(:))'; 
        nb=size(b); nb=nb(2);
        b=[zeros(1,(d-nb)) b]; 
    c=dec2bin(k);
        c=str2num(c(:))'; 
        nc=size(c); nc=nc(2); 
        c=[zeros(1,(d-nc)) c]; 
    eps=0.0000001;
    if(pf(a,b)<0.5-eps)&&(pf(a,c)<0.5-eps)&&(pf3(a,b,c)>pf3(b,a,c))&&(pf3(a,b,c)>pf3(c,a,b))
        fprintf('%2d %2d %2d %f %f %f | %f %f %f', ...
            i,j,k, pf(a,b), pf(a,c), pf(b,c), pf3(a,b,c), pf3(b,a,c),pf3(c,a,b));
        if ((pf(b,c)>0.5+eps)&& (pf3(c,a,b)>pf3(b,a,c)))
            fprintf('*');
        end
        if ((pf(b,c)<0.5-eps)&& (pf3(c,a,b)<pf3(b,a,c)))
            fprintf('*');
        end
        fprintf('\n');
    end
    
end
end
end
end

function y=exp0(a)

global p; 
global q;
 
n=size(a);
n=n(2);

exp=0; 
exp=exp+p^(-sum(a))*q^(-sum(1-a));
for (i=2:n)
    a1=a(i:n);
    a2=a(1:(n-i+1));
    if (isequal(a1,a2))
        exp=exp+p^(-sum(a1))*q^(-sum(1-a1));
    end
end
y=exp;
end

function y=exp1(a,b)

global p; 
global q;
 

na=size(a);
na=na(2);

nb=size(b);
nb=nb(2);

done=(na==0);
i=1;
exp=exp0(b);
while (done==0)
    if (na-i+1<nb)
        a2=a(i:na);
        b1=b(1:(na-i+1));
        if (isequal(b1,a2))
            exp=exp-exp0(b1);
            done=1; 
        end
    end
    i=i+1;
    if (i>na) done=1;
    end
end

y=exp;
 
end

function y=pf(a,b)

M=[0 1 1; 1 0 exp1(b,a); 1 exp1(a,b) 0];
R=[1 exp0(a) exp0(b)]';
vp=M^(-1)*R; 
y=vp(2);

end

function y=pf3(a,b,c)

M=[0 1 1 1; 1 0 exp1(b,a) exp1(c,a); 1 exp1(a,b) 0 exp1(c,b); 1 exp1(a,c) exp1(b,c) 0];
R=[1 exp0(a) exp0(b) exp0(c)]';
vp=M^(-1)*R; 
y=vp(2);

end
