Linear problem solving with matrix constraints in Rust

Question

I am trying to rewrite a fairness ranking algorithm (source: https://arxiv.org/abs/1802.07281) from Python to Rust. The objective is finding a document-ranking probability matrix that is doubly stochastic and, by use of an utility vector (i.e. the document relevance in this case) gives fair exposure to all document types.

The objective is thus to maximise the expected utility under the following constraints:

• sum of probabilities for each position equals 1;
• sum of probabilities for each document equals 1;
• every probibility is valid (i.e. 0 <= P[i,j] <= 1);
• P is fair (disparate treatment constraints).

In Python we have done this using CVXPY:

    u = documents[['rel']].iloc[:n].values.ravel()  # utility vector
    v = np.array([1.0 / (np.log(2 + i)) for i in range(n)])  # position discount vector
    P = cp.Variable((n, n))  # linear maximizing problem uͭPv s.t. P is doubly stochastic and fair.
    
    # Construct f in fͭPv such that for P every group's exposure divided by mean utility should be
    # equal (i.e. enforcing DTC). Do this for the set of every individual two groups:
    #                 example: calculated f for three groups {a, b, c}
    #   resulting constraints: [a - b == 0, a - c == 0, b - c == 0]
    groups = {k: group.index.values for k, group in documents.iloc[:n].groupby('document_type')}
    fairness_constraints = []
    for k0, k1 in combinations(groups, 2):
        g0, g1 = groups[k0], groups[k1]
        f_i = np.zeros(n)
    
        f_i[g0] = 1 / u[g0].sum()
        f_i[g1] = -1 / u[g1].sum()
        fairness_constraints.append(f_i)
    
    # Create convex problem to solve for finding the probabilities that
    # a document is at a certain position/rank, matching the fairness criteria
    objective = cp.Maximize(cp.matmul(cp.matmul(u, P), v))
    constraints = ([cp.matmul(np.ones((1, n)), P) == np.ones((1, n)),  # ┤
                    cp.matmul(P, np.ones((n,))) == np.ones((n,)),      # ┤
                    0.0 <= P, P <= 1] +                                # └┤ doubly stochastic matrix constraints
                   [cp.matmul(cp.matmul(c, P), v) == 0 for c in fairness_constraints])  # DTC
    
    prob = cp.Problem(objective, constraints)
    prob.solve(solver=cp.CBC)

This works great for multiple solvers, including SCS, ECOS and CBC.

Now trying to implement the algorithm above to Rust, I have resolved to crates like good_lp and lp_modeler. These should both be able to solve linear problems using CBC as also demonstrated in the Python example above. I am struggling however to find examples on how to define the needed constraints on my matrix variable P.

The code below is my work in progress for rewriting the Python code in Rust, using in this case the lp_modeler crate as an example. The code below compiles but panics when run. Furthermore I don't know how to add the disparate treatment constraints in a way Rust likes, as no package seems to be able to accept equality constraints on two vectors.

let n = cmp::min(u.len(), 25);

let u: Array<f32, Ix1> = array![...];  // utility vector filled with dummy data

// position discount vector
let v: Array<f32, Ix1> = (0..n)
    .map(|i| 1.0 / ((2 + i) as f32).ln())
    .collect();

let P: Array<f32, Ix2> = Array::ones((n, n));

// dummy data for document indices and their types
let groups = vec![
    vec![23], // type A
    vec![8, 10, 16, 19], // type B
    vec![0, 1, 2, 3, 4, 5, 6, 7, 9, 11, 12, 13, 15, 21, 24], // type C
    vec![14, 17, 18, 20, 22] // type D
];

let mut fairness_contraints: Vec<Vec<f32>> = Vec::new();
for combo in groups.iter().combinations(2).unique() {
    let mut f_i: Vec<f32> = vec![0f32; n];

    {  // f_i[g0] = 1 / u[g0].sum()
        let usum_g0: f32 = combo[0].iter()
            .map(|&i| u[i])
            .sum();

        for &i in combo[0].iter() {
            f_i[i] = 1f32 / usum_g0;
        }
    }

    {  // f_i[g1] = -1 / u[g1].sum()
        let usum_g1: f32 = combo[1].iter()
            .map(|&i| u[i])
            .sum();

        for &i in combo[1].iter() {
            f_i[i] = -1.0 / usum_g1;
        }
    }

    fairness_contraints.push(f_i);
}

let mut problem = LpProblem::new("Fairness", LpObjective::Maximize);

problem +=  u.dot(&P).dot(&v); // Expected utility objective

// Doubly stochastic constraints
for col in P.columns() {  // Sum of probabilities for each position
    problem += sum(&col.to_vec(), |&el| el).equal(1);
}

for row in P.rows() {  // Sum of probabilities for each document
    problem += sum(&row.to_vec(), |&el| el).equal(1);
}

// Valid probability constraints
for el in P.iter() {
    problem += lp_sum(&vec![el]).ge(0);
    problem += lp_sum(&vec![el]).le(1);
}

// TODO: implement DTC fairness constraints

let solver = CbcSolver::new();
let result = solver.run(&problem);

Can anybody give me a nudge in the right direction on this specific problem? Thanks in advance!