function  [v, S, r, ri, mi] = cbfq(q, p, d)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% File:               cbfq.m
%
% Author:             Tasos Kalandrakis
%
% Description:        Computes equilibrium values of Baron and Ferejohn
%                     (1989) model with general quota voting rule,
%                     recognition probabilities and discount factors.
%
% Created:            Feb - 01 - 2013
%
% Last Modified:      Oct - 16 - 2022
%
% Language:           MATLAB
%
% Related Reference:  "Computation of equilibrium values in the Baron and
%                     Ferejohn bargaining model," Games and Economic
%                     Behavior, forthcoming.
%
% Acknowledgments:    Thanks to Jan Zapal and to Ognyan Kounchev for 
%                     alerting me to a couple of bugs; any remaining errors
%                     are mine.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% INPUT
% q is quota (1<= q<= n)
% p is 1 x n vector of recognition probabilities
% d is 1 x n vector of discount factors

% OUTPUT
% v is 1 x n vector of equilibrium values
% S is maximum surplus
% r is marginal reservation value
% ri is 1 x n vector of reservation values
% mi is 1 x n vector of probabilities of inclusion in winning coalition


%check recognition probabilities
if sum(p)~=1
    fprintf('\n')
    fprintf('WARNING: recognition probabilities do not sum up to 1. Applying proportional adjustment.\n')
    p=p/sum(p);
end

%preliminaries
dp=d.*p;
Lr=dp./(1-d);
Hr=dp./(1-dp);
Lv=(p./(1-d));

overTh=sort(Lr);
overTh=overTh(q);

underTh=sort(Hr);
underTh=underTh(q);

Theta=union(Lr,Hr);
Theta=sort(Theta(((Theta>=underTh)&(Theta<=overTh))),'descend');
Theta=[Theta(2:end)' Theta(1:(end-1))'];
K=length(Theta(:,1));

invd=1./d;
invd((d==0))=0;
b=ones(2,1);
Theta(:,3)=(Theta(:,1)+Theta(:,2))./2;


%main algorithm
if overTh==underTh
    h=bsxfun(@le,overTh,Hr);
    l=bsxfun(@gt,overTh,Lr);
    m=bsxfun(@gt,overTh,Hr) & bsxfun(@le,overTh,Lr);
    S=(1+sum(l.*Lr))/(1+overTh*(q-sum(l)));
    r=overTh * S;
else
    %initial value
    Kset=zeros(K,1);
    ratio=Theta(ceil(K/2),3);
    h=bsxfun(@le,ratio,Hr);
    l=bsxfun(@gt,ratio,Lr);
    m=bsxfun(@gt,ratio,Hr) & bsxfun(@le,ratio,Lr);
    A(1,1)=1+sum(l.*Lr);
    A(1,2)=q-sum(l);
    A(2,1)=sum(l.*Lv)+sum(h.*p);
    A(2,2)=sum(m.*invd)+sum(h.*p);
    x=A\b;
    S=x(1);
    r=x(2);
    k=(zeros(K,1)==1);
    k(ceil(K/2),1)=(1==1);
    
    fprintf('\n')
    while ((r<Theta(k,1)*S)||(r>Theta(k,2)*S))
        fprintf('\n')
        if ((r-Theta(k,1)*S>-eps)&&(r-Theta(k,2)*S<eps))
            fprintf('WARNING: iterate outside piece of linearity likely due to machine precision.\n')
            break
        end            
        Kset=Kset|k;
        check=bsxfun(@ge,r,Theta(:,1).*S)&bsxfun(@le,r,Theta(:,2).*S)&(~Kset);
        if sum(check)==0
            ang=Theta(~Kset,3);
            [~,I]=min(abs(ang-ratio));
            ratio=ang(I);
            if ((r>overTh*S)||(r<underTh*S))
                fprintf('Iteration %u. Iterate out of bounds. Restart at nearest piece of linearity.\n', sum(Kset))
            else
                fprintf('Iteration %u. A cycle was detected. Restart at nearest piece of linearity.\n', sum(Kset))
            end
        else
            ang=Theta(check,3);
            ratio=min(ang);
            fprintf('Iteration %u.\n',sum(Kset))
        end
        h=bsxfun(@le,ratio,Hr);
        l=bsxfun(@gt,ratio,Lr);
        m=bsxfun(@gt,ratio,Hr) & bsxfun(@le,ratio,Lr);
        A(1,1)=1+sum(l.*Lr);
        A(1,2)=q-sum(l);
        A(2,1)=sum(l.*Lv)+sum(h.*p);
        A(2,2)=sum(m.*invd)+sum(h.*p);
        res=norm(A*[S; r]-b);
        fprintf('     Residual: %g\n',res)
        x=A\b;
        S=x(1);
        r=x(2);
        k=bsxfun(@gt,ratio,Theta(:,1))&bsxfun(@lt,ratio,Theta(:,2));
    end
    res=norm(A*[S; r]-b);
    
    fprintf('\n')
    fprintf('Convergence achieved in %u iterations.  Residual is: %g (may not be zero due to machine precision)\n',sum(Kset),res)
end

%Final computations
v=l.*(Lv.*S)+(m.*invd)*r +(h.*p)*(S+r);
ri=l.*(Lr.*S)+m.*r +h.*(Hr.*S);
if r>0
    mi=l.*(1-p)+m.*(invd-p.*(S/r)-p);
else
    mi=l.*(1-p);
end

end