Continuous knapsack problem in Rust

Question

Background

The continuous knapsack problem is the following linear program:

$$ \begin{align} \text{maximize} \quad & f(x) = \sum_{i=1}^n u_i x_i \\ \text{subject to} \quad & \sum_{i=1}^n w_i x_i \leq W \\ & 0 \leq x_i \leq 1, ~~ i = 1 \dots n \end{align} $$

As discussed in the link above, it is solvable in \$O(n\log n)\$ time. First you index the items in descending order by \$u_i / w_i\$ ^{(edit: I got this backwards in my implementation but will leave it as-is to preserve the correctness of preexisting answers)}, then:

For each material, the amount \$x_i\$ is chosen to be as large as possible:

• If the sum of the choices made so far equals the capacity \$W\$, then the algorithm sets \$x_i = 0\$.
• If the difference \$d\$ between the sum of the choices made so far and \$W\$ is smaller than \$w_i\$, then the algorithm sets \$x_i = d / w_i\$.
• In the remaining case, the algorithm chooses \$x_i = 1\$.

Motivation

I primarily use Julia and Python in my work, but I am trying to learn Rust and C++ because many of the tasks I need to do are simple numerical schemes like this that run much faster when coded in a lower-level language.

Eventually, I am hoping to try implementing a branch-and-bound scheme for integer programming in Rust, but for now I decided to start with this problem.

Questions

The first wrinkle I encountered here was the lack of an argsort function in Rust. In the link above, it is suggested to use the sort_by_key() function, but this doesn't work for floats because they are not totally ordered. I opted instead to create an index set and sort_by() it in place using a closure that references the array I wanted to argsort on. Is this the best approach?

I also sense that my choice to use arrays for the input data and vectors for the numerical manipulation was somewhat arbitrary. In principle, once we have read in the data, we know its size, so it should be possible to store r and x below as n-arrays instead. However, trying to initialize let mut x = [0.0; n] throws an error because n is not a constant, and I was unable to wrangle it into one. Is there a way to use arrays instead of vectors below? Would it be wise to do so?

Finally, I have a somewhat nebulous sense of how I should handle the I/O for an algorithm like this in an industrial context. In Julia, for example, my workflow would be to put my solution algorithms in a module, then put my problem data in a CSV in the same directory as the module, and run julia --project in that directory so I could use another interface like DataFrames.jl to pass the data from the CSV to my function.

In Rust, there is no REPL (obviously), so it seems like we are stuck with either hardcoding the problem data into the main() function as I have done below (which requires recompiling the code) or putting some read_line()s in the main() function (which is cumbersome). Are there any examples of Rust programs that work somewhat like the Julia idea described above—compile the code once, then solve arbitrary problems saved as a CSV or text file in a local directory?

Finally, a general code review would be more than welcome.

Implementation

Code:

fn main() {
    let u = [3, 6, 4, 7, 8];    // Utility values
    let w = [1, 4, 2, 5, 7];    // Weights
    let w_max = 14;             // Knapsack capacity

    let (x, f) = solve_continuous_knapsack(&u, &w, w_max);

    println!("  x  = {:?}", x);
    println!("f(x) = {}", f);
}

fn solve_continuous_knapsack(u: &[isize], w: &[isize], w_max: isize) -> (Vec<f64>,f64) {
    // Check that the problem is well-posed
    let n = u.len();
    assert_eq!(n, w.len());
    for i in 0..n {assert!(w[i] > 0)};

    // Utility/cost ratio for each item
    let mut r: Vec<f64> = Vec::new();

    // Container for continuous solution
    let mut x: Vec<f64> = Vec::new();

    for i in 0..n {
        r.push(u[i] as f64 / w[i] as f64);
        x.push(0.0);
    }

    // Sort the items in ascending order by r[i]
    let mut index_set = (0..n).collect::<Vec<usize>>();
    index_set.sort_by(|&i, &j| (&r[i]).partial_cmp(&r[j]).unwrap());

    // Iterate through the indices until the knapsack is full
    let mut f = 0.0;
    let mut w_left = w_max;
    for &i in index_set.iter() {
        if w[i] <= w_left {
            x[i] = 1.0;
            f += u[i] as f64;
            w_left -= w[i];
        } else {
            x[i] = w_left as f64 / w[i] as f64;
            f += u[i] as f64 * x[i];
            break;
        }
    }

    return (x, f);
}

Output:

  x  = [0.0, 0.5, 0.0, 1.0, 1.0]
f(x) = 18

Answer 1

Is this the best approach?

My approach here is to use the noisy_float crate. It provides types that wrap floats, asserting that they do not contain NaN and being totally ordered.

Is there a way to use arrays instead of vectors below? Would it be wise to do so?

You can only use arrays when the size is known at compile time. In any other case, just use a Vec. Vec does not have any extra overhead compared to a runtime sized array in other languages, so there's no reason to want that in Rust.

Are there any examples of Rust programs that work somewhat like the Julia idea described above—compile the code once, then solve arbitrary problems saved as a CSV or text file in a local directory?

In such cases, my preference is to take input from a JSON file which I read in via the serde_json crate. This way I can simply define the input format as a struct and read it and process it really easily. If I can't control the format, I typically find a crate that knows how to read the format. For example there is a csv create that parses csv input.

fn solve_continuous_knapsack(u: &[isize], w: &[isize], w_max: isize) -> 
(Vec<f64>,f64) {

I recommend taking the inputs as f64 instead of isize, so you convert to f64 in (almost?) all cases anyways.

// Utility/cost ratio for each item
let mut r: Vec<f64> = Vec::new();

// Container for continuous solution
let mut x: Vec<f64> = Vec::new();

for i in 0..n {
    r.push(u[i] as f64 / w[i] as f64);
    x.push(0.0);
}

This strategy of creating an empty Vec and then pushing into won't be the most performant. You're better off collecting from an iterator.

let mut r: Vec<f64> = u.iter().zip(w.iter()).map(|(&u, &w)| (u as f64) / (w as f64)).collect();
let mut x: Vec<f64> = u.iter().map(|_| 0.0).collect();