function [R r]=triangle_ob_lp(d1,d1E,d2,d2E,d3,d3E)
%[R r]=triangle_ob_lp(d1,d1E,d2,d2E,d3,d3E) calculates the outer bound
%for the triangle network from the channel erasure probabilities. 
%
%
%The input parameters are standing for the erasure 
%probabilities:
% - of the S-D channel (d1), 
% - of Eve on the S-D channel (d1E), 
% - of the S-U channel (d2), 
% - of Eve on the S-U channel (d2E),
% - of the U-D channel (d3),
% - of Eve on the U-D channel (d3E)
%in this order.
%
%Output R is the value of the outer bound LP, while r is the solution vector.
%The index-variable map of r is available in a separate document.
%

%Number of variables
v = 41;

%Objective function. The last variable is the rate. Since Matlab solves a minimazation problem the coefficent here is -1
%meaning that we seek the largest rate.
f = zeros(v,1);
f(v) = -1;

%We collect constraints in matrix A. Vector b collects the constants of the inequalities (mostly 0). Aeq and beq correspond to equality constraints.
b = zeros(0,1);
A = zeros(0,v);
Aeq = zeros(0,v);

%Rate constraints from the S-UD cut
%inequalities (101)-(103) from the proof
A = [A;zeros(1,v)]; A(end,[v 1 2 ]) = [1 -(1-d1) -(1-d2) ];
A = [A;zeros(1,v)]; A(end,[v 3 4 ]) = [1 -(1-d1*d1E) -(1-d2)  ];
A = [A;zeros(1,v)]; A(end,[v 5 6 ]) = [1 -(1-d2*d2E) -(1-d1) ];

%Rate constraints from the SU-D cut
%inequalities (104)-(106)
A = [A;zeros(1,v)]; A(end,[v 7 8  ]) = [1 -(1-d1) -(1-d3) ];
A = [A;zeros(1,v)]; A(end,[v 9 10  ]) = [1 -(1-d1*d1E) -(1-d3)  ];
A = [A;zeros(1,v)]; A(end,[v 11 12  ]) = [1 -(1-d3*d3E) -(1-d1) ];

%Key generation constraints 
%inequalities (107)-(110)
A = [A;zeros(1,v)]; A(end,[13 14 15 ]) = [(1-d1E) -d1E*(1-d1) -(1-d2) ];
A = [A;zeros(1,v)]; A(end,[16 17 18 ]) = [(1-d2E) -d2E*(1-d2) -(1-d1) ];

A = [A;zeros(1,v)]; A(end,[19 20 21 ]) = [(1-d3E) -d3E*(1-d3) -(1-d1) ];
A = [A;zeros(1,v)]; A(end,[22 23 24 ]) = [(1-d1E) -d1E*(1-d1) -(1-d3) ];

%Security on SU-D cut
%inequalities (111)-(112)
A = [A;zeros(1,v)]; A(end,[3 13]) = [1 -1];
A = [A;zeros(1,v)]; A(end,[5 16]) = [1 -1];

%Security on SD-U cut
%inequalities (113)-(114)
A = [A;zeros(1,v)]; A(end,[11 19]) = [1 -1];
A = [A;zeros(1,v)]; A(end,[9 22]) = [1 -1];

%Key reuse inequalities
%inequalities (115)-(116)
A = [A;zeros(1,v)]; A(end,[25 23 26]) = [(1-d1E) -d1E*(1-d1) -d1E*(1-d1)];
A = [A;zeros(1,v)]; A(end,[27 20 28]) = [(1-d3E) -d3E*(1-d3) -d3E*(1-d3)];

%Equalities (117)-(118)
Aeq = [Aeq;zeros(1,v)]; Aeq(end,[22 25 26]) = [-1 1 1];
Aeq = [Aeq;zeros(1,v)]; Aeq(end,[19 27 28]) = [-1 1 1];

%%Connecting the two cuts
%inequalities (119)-(123)
A = [A;zeros(1,v)]; A(end,[15 29 30 ]) = [-(1-d2) (1-d1*d1E) (1-d3*d3E) ];
A = [A;zeros(1,v)]; A(end,[28 31 32 33 29  ]) = [(1-d3*d3E) (1-d3*d3E) -(1-d2) -(1-d2) -(1-d1*d1E)  ];
A = [A;zeros(1,v)]; A(end,[26 30 34 35 ]) = [(1-d1*d1E) -(1-d3) -(1-d3) (1-d3) ];

A = [A;zeros(1,v)]; A(end,[36 34 35 30 ]) = [(1-d1*d1E) (1-d3) -(1-d3) -d3*(1-d3E) ];
A = [A;zeros(1,v)]; A(end,[32 14 37 ]) = [(1-d2) -(1-d1)  -(1-d1) ]; 
%equality (124)
Aeq = [Aeq;zeros(1,v)]; Aeq(end,[14 38 29]) = [-1 1 -1];

%Trivial equalities
%equalities (125)-(128)
Aeq = [Aeq;zeros(1,v)]; Aeq(end,[23 38 36]) = [-1 1 1];
Aeq = [Aeq;zeros(1,v)]; Aeq(end,[20 39 31]) = [-1 1 1];
Aeq = [Aeq;zeros(1,v)]; Aeq(end,[21 38 40]) = [-1 1 1];
Aeq = [Aeq;zeros(1,v)]; Aeq(end,[24 39 34]) = [-1 1 1];

%Transmissions
%inequalities (130)-(141)
A = [A;zeros(1,v)]; A(end,[18 6]) = [1 1];
b(size(A,1)) = 1;
A = [A;zeros(1,v)]; A(end,[15 4]) = [1 1];
b(size(A,1)) = 1;
A = [A;zeros(1,v)]; A(end,[14 1 37]) = [1 1 1];
b(size(A,1)) = 1;
A = [A;zeros(1,v)]; A(end,[15 2 33 ]) = [1 1 1];
b(size(A,1)) = 1;
A = [A;zeros(1,v)]; A(end,[18 1]) = [1 1 ];
b(size(A,1)) = 1;
A = [A;zeros(1,v)]; A(end,[17 2 ]) = [1 1 ];
b(size(A,1)) = 1;
A = [A;zeros(1,v)]; A(end,[21 12]) = [1 1];
b(size(A,1)) = 1;
A = [A;zeros(1,v)]; A(end,[24 10]) = [1 1];
b(size(A,1)) = 1;
A = [A;zeros(1,v)]; A(end,[23 7]) = [1 1];
b(size(A,1)) = 1;
A = [A;zeros(1,v)]; A(end,[24 8]) = [1 1];
b(size(A,1)) = 1;
A = [A;zeros(1,v)]; A(end,[21 7]) = [1 1];
b(size(A,1)) = 1;
A = [A;zeros(1,v)]; A(end,[20 8]) = [1 1];
b(size(A,1)) = 1;

%Pad with zeros if no constants were given
if length(b) < size(A,1)
    b(size(A,1)) = 0;
end

%Every variable is nonnegative
A = [A; -eye(v)];
b = [b zeros(1,v)];

%No equality has any constant
beq = zeros(size(Aeq,1),1);

r = linprog(f,A,b',Aeq,beq);
R = r(v);