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