(*Code to generate input for the Fourier-Motzkin solver PORTA in order to generate all Shannon-type entropic inequalities in a marginal scenario*) (*Equivalently code to generate ncessary and sufficient conditions for non-contextuality of a partial polymatroid*)
(*Written by Rafael Chaves (rafael.chaves@icfo.es) and Tobias Fritz (tobias.fritz@icfo.es) in 2011*)

(*The user needs to specify the number "n" of random variables {(Subscript[A, 1]),...,(Subscript[A, n])} and also the marginal scenario defined by the list "M"*)

(*In this example we consider the marginal scenario for the n-cycle described in the paper with n=4*)
n=4;

(*For the marginal scenario, one only needs to specify the maximal subsets*)
M={{1,2},{2,3},{3,4},{4,1}};

(*A PORTA variable "xi" for 1≤ i ≤(2^n) represents the coefficient of a certain joint entropy H((Subscript[A, S])). The number "i" in binary representation specifies the characteristic function of the set of variables S.*)
(*For example "x4" corresponds to H((Subscript[A, n-1])(Subscript[A, n])). The following function "getPortaStr" generates the list of those PORTA variables which should be eliminated by the Fourier-Motzkin procedure. These correspond to those joint entropies H((Subscript[A, S])) for S not in M.*)

getPortaStr[M_,n_]:=Module[{csets,k,i,s,res},
(*initialise*)
res={};
i=0;
(*all compatible subsets*)
csets=Flatten[Subsets[#,{1,n}]&/@M,1];

(*loop through positions in output string*)
For[k=1,k≤(2^n),k++,
(*convert position to integer form*)
s=IntegerDigits[k-1,2,n];
s=Flatten[Position[s,1]];
(*check if compatible subset exists*)
AppendTo[res,If[MemberQ[csets,s],0,i++;i]];
];

StringJoin[Riffle[ToString/@res," "]]
]

(*Start the string that will be exported as a .ieq file (input for PORTA)*)
str= {"DIM="<>ToString[(2^n)],"\n","ELIMINATION_ORDER",getPortaStr[M,n],"\n"};

str=Join[str,{"INEQUALITIES_SECTION"}];

(*The following loop generates all the submodularity inequalities of Eq. (8) in the manuscript*)
For[i=(2^n)-1,i≥0,i--,
For[j=i-1,j≥0,j--,
aux1=0;
aux2=0;
aux3=0;
aux4=0;
aux5=0;

(*Calculate the Hamming distance to assure that H(S U i) and H(S U j) differ in at most one position*)
aux5=HammingDistance[IntegerDigits[i,2,n],IntegerDigits[j,2,n]];

(*Calculate if H(S U j) == H(S) or not *)
If[j== BitAnd[i,j],aux1=1,aux1=0];

(*Calculate if H(S U j) == H(S U {i,j}) or not *)
If[j== BitOr[i,j],aux2=1,aux2=0];

(*Calculate if H(S U i) == H(S) or not *)
If[i== BitAnd[i,j],aux3=1,aux3=0];

(*Calculate if H(S U i) == H(S U {i,j}) or not *)
If[i== BitOr[i,j],aux4=1,aux4=0];

(*If H(S U i) and H(S U j) differ in at most one position, the submodularity inequality is included*)
If[aux5≤2,If[BitOr[BitOr[BitOr[aux1,aux2],aux3],aux4]==0,

AppendTo[str,"x"<>ToString[i+1]<>" + "<>"x"<>ToString[j+1]<>" >=  x"<>ToString[BitOr[i,j]+1]<>" + x"<>ToString[BitAnd[i,j]+1]]; 

]]]]


(*The following loop generates all the inequalities f([n]\{i}) ≤ f([n]) for all i ϵ [n]*)
vec10=Permutations[Table[If[i≠ n,1,0],{i,n}]];
For[i=1,i≤n,i++,AppendTo[str,"x"<>ToString[(2^n)]<>" >= x"<>ToString[FromDigits[vec10[[i]],2]+1]]]

(*Generate the last constraint, namely that f(∅)=0 (∅=empty set)*)
str=Join[str,{"x1 == 0"}];

(*Any additional constraint, for example conditional independence relations, should be amended here by hand.*)
(*This example includes the independence between (Subscript[A, n-1]) and (Subscript[A, n]) as an aditional constraint that should be written as "x2 + x3 - x4 == 0", corresponding to H((Subscript[A, n-1]))+H((Subscript[A, n]))=H((Subscript[A, n-1])(Subscript[A, n])).*)
str=Join[str,{"x2 + x3 - x4 == 0"}];

(*This ends the file*)
str=Join[str,{"\n","END","\n"}]

(*Export the string to a text file in the format used by PORTA. Here the user need to define a path and a name file. Then run the PORTA command "fmel" on the file.*)
Export["C:\\inputporta.ieq",str,"Text"]