import numpy as np
from scipy.optimize import linprog

#potentials
f = np.array([1., 0., 0., 1., 0., 5., 0., 0.])

#local probability constraints
##############################
## specify the constraints which have a value of one as their right hand side
## Dimensions: A_eq1 (list of three lists each with 8 entries)
##############################
A_eq1 = [[1,1 , 0, 0, 0, 0, 0,0 ],[0,0 ,1 ,1 ,0 ,0 ,0 ,0 ],[0,0 ,0 ,0 ,1 ,1 ,1 ,1 ]]
b_eq1 = [[1],]*3

#marginalization constraints
##############################
## specify the constraints which have a value of zero as their right hand side
## Dimensions: A_eq2 (list of three lists each with 8 entries)
##############################
A_eq2 = [[-1, 0, 0,0 ,1 ,0 , 1, 0],[0, -1, 0, 0, 0, 1,0 , 1],[0, 0, -1, 0, 1, 1, 0,0 ],[0, 0, 0, -1,0 ,0 ,1 ,1 ]]
b_eq2 = [[0],]*4

#bounds
bounds = [(0, 1),]*8

res = linprog(-f, A_eq=np.concatenate((A_eq1,A_eq2)), b_eq=np.concatenate((b_eq1,b_eq2)), bounds=bounds)
print(res)

print(np.matmul(A_eq1, res.x))
print(np.matmul(A_eq2, res.x))
043d708eb533dde1f9a724a079fd4ff0c3c243ac
<mask token>

<mask token>
print(res)
print(np.matmul(A_eq1, res.x))
print(np.matmul(A_eq2, res.x))

<mask token>
f = np.array([1.0, 0.0, 0.0, 1.0, 0.0, 5.0, 0.0, 0.0])
A_eq1 = [[1, 1, 0, 0, 0, 0, 0, 0], [0, 0, 1, 1, 0, 0, 0, 0], [0, 0, 0, 0, 1,
    1, 1, 1]]
b_eq1 = [[1]] * 3
A_eq2 = [[-1, 0, 0, 0, 1, 0, 1, 0], [0, -1, 0, 0, 0, 1, 0, 1], [0, 0, -1, 0,
    1, 1, 0, 0], [0, 0, 0, -1, 0, 0, 1, 1]]
b_eq2 = [[0]] * 4
bounds = [(0, 1)] * 8
res = linprog(-f, A_eq=np.concatenate((A_eq1, A_eq2)), b_eq=np.concatenate(
    (b_eq1, b_eq2)), bounds=bounds)
print(res)
print(np.matmul(A_eq1, res.x))
print(np.matmul(A_eq2, res.x))

import numpy as np
from scipy.optimize import linprog
f = np.array([1.0, 0.0, 0.0, 1.0, 0.0, 5.0, 0.0, 0.0])
A_eq1 = [[1, 1, 0, 0, 0, 0, 0, 0], [0, 0, 1, 1, 0, 0, 0, 0], [0, 0, 0, 0, 1,
    1, 1, 1]]
b_eq1 = [[1]] * 3
A_eq2 = [[-1, 0, 0, 0, 1, 0, 1, 0], [0, -1, 0, 0, 0, 1, 0, 1], [0, 0, -1, 0,
    1, 1, 0, 0], [0, 0, 0, -1, 0, 0, 1, 1]]
b_eq2 = [[0]] * 4
bounds = [(0, 1)] * 8
res = linprog(-f, A_eq=np.concatenate((A_eq1, A_eq2)), b_eq=np.concatenate(
    (b_eq1, b_eq2)), bounds=bounds)
print(res)
print(np.matmul(A_eq1, res.x))
print(np.matmul(A_eq2, res.x))

import numpy as np
from scipy.optimize import linprog

#potentials
f = np.array([1., 0., 0., 1., 0., 5., 0., 0.])

#local probability constraints
##############################
## specify the constraints which have a value of one as their right hand side
## Dimensions: A_eq1 (list of three lists each with 8 entries)
##############################
A_eq1 = [[1,1 , 0, 0, 0, 0, 0,0 ],[0,0 ,1 ,1 ,0 ,0 ,0 ,0 ],[0,0 ,0 ,0 ,1 ,1 ,1 ,1 ]]
b_eq1 = [[1],]*3

#marginalization constraints
##############################
## specify the constraints which have a value of zero as their right hand side
## Dimensions: A_eq2 (list of three lists each with 8 entries)
##############################
A_eq2 = [[-1, 0, 0,0 ,1 ,0 , 1, 0],[0, -1, 0, 0, 0, 1,0 , 1],[0, 0, -1, 0, 1, 1, 0,0 ],[0, 0, 0, -1,0 ,0 ,1 ,1 ]]
b_eq2 = [[0],]*4

#bounds
bounds = [(0, 1),]*8

res = linprog(-f, A_eq=np.concatenate((A_eq1,A_eq2)), b_eq=np.concatenate((b_eq1,b_eq2)), bounds=bounds)
print(res)

print(np.matmul(A_eq1, res.x))
print(np.matmul(A_eq2, res.x))