Borrowing problem from porting some C++ code

help
1

I have a boring common borrowing problem, and I'd like to know if you have suggestions and improvements.

I am trying to translate to Rust this little C++ program:
https://sourceforge.net/p/ctidbits/code/HEAD/tree/trunk/discrepancy/discrepancy.cpp

The input file is:
https://sourceforge.net/p/ctidbits/code/HEAD/tree/trunk/discrepancy/discrepancy.dat?force=True

The C++ code (with tiny changes):

#include <fstream>
#include <vector>
#include <algorithm>
#include <cmath>

struct Point {
    double x, y;
    bool marked;
    Point(double x_, double y_) : x(x_), y(y_), marked(false) {}
    Point(const Point &p) : x(p.x), y(p.y), marked(p.marked) {}
};

std::vector<Point> points;
std::vector<Point*> sortX, sortY, subSetY;

double maxDiscrepancy, recMinX, recMinY, recMaxX, recMaxY;

void readPoints() {
    std::ifstream inp("discrepancy.dat");
    double x, y;
    while (inp >> x >> y)
        points.push_back(Point(x, y));
    inp.close();
    printf("Read %u points from input file\n", points.size());
}

bool pointCompX(const Point *a, const Point *b) { return a->x < b->x; }
bool pointCompY(const Point *a, const Point *b) { return a->y < b->y; }

void sortPoints() {
    sortX.resize(points.size());
    sortY.resize(points.size());
    for (size_t i = 0; i < points.size(); i++) {
        sortX[i] = &points[i];
        sortY[i] = &points[i];
    }

    std::sort(sortX.begin(), sortX.end(), pointCompX);
    std::sort(sortY.begin(), sortY.end(), pointCompY);
}

void innerLoopY(int n, double lowX, double highX) {
    // printf("inner loop %d\n", n);
    const double xDim = highX - lowX;
    for (int i = 0; i < n; ++i) {
        double lowY = subSetY[i]->y;
        for (int j = i; j < n; ++j) {
            const double highY = subSetY[j]->y;
            const double yDim = highY - lowY;
            const int count = j + 1 - i;
            const double Discrepancy = fabs(points.size() * xDim * yDim - count);

            if (Discrepancy > maxDiscrepancy) {
                maxDiscrepancy = Discrepancy;
                recMinX = lowX;
                recMinY = lowY;
                recMaxX = highX;
                recMaxY = highY;
            }
        }
    }
}

void outerLoopX() {
    subSetY.resize(points.size());
    maxDiscrepancy = 0.0;

    for (size_t i = 0; i < points.size(); ++i) {
        printf("Low X = %u\n", i);
        double lowX = sortX[i]->x;
        for (size_t j = i; j < points.size(); ++j) {
            // printf("High X = %u\n", j);
            sortX[j]->marked = true;
            double highX = sortX[j]->x;
            int l = 0;
            for (size_t k = 0; k<points.size(); ++k)
                if (sortY[k]->marked) {
                    subSetY[l] = sortY[k];
                    l++;
                }
            innerLoopY(l, lowX, highX);
        }

        for (size_t j = i; j < points.size(); ++j)
            sortX[j]->marked = false;
    }
}

int main() {
    readPoints();
    sortPoints();
    outerLoopX();

    printf("Maximum Discrepancy = %f\n", maxDiscrepancy);
    printf("Rect %f %f %f %f\n", recMinX, recMinY, recMaxX, recMaxY);
}

As you see it contains four arrays (vectors): pointlist that contains the input data points. Its x and y coordinates never change once loaded, only Point.marked changes. There are also sortX and sortY arrays that contain pointers to pointlist (and this causes borrowing problems in Rust) that point to Points that are sorted according to the x or y field. And there's subSetY, another array of pointers that are sorted by x field, and are a shrinking selection (subset) of the Points.

Minutes after starting a Rust translation of the C++ code I've seen the borrowing problems, so I've soon modified the algorithm, replacing the pointers with array indexes, and I have also removed the global variables to make the Rust code nicer (equally removing the global variables in the C++ code makes it a little faster):

struct Point {
    x: f64,
    y: f64,
    marked: bool
}

fn read_data() -> Vec<Point> {
    use std::fs::File;
    use std::io::Read;

    let mut data = String::new();
    File::open("discrepancy.dat").unwrap().read_to_string(&mut data).unwrap();

    let points: Vec<Point> =
        data
        .lines()
        .map(|line| line.split(' ').map(|p| p.parse().unwrap()))
        .map(|mut parts| Point { x: parts.next().unwrap(),
                                 y: parts.next().unwrap(),
                                 marked: false })
        .collect();

    println!("Read {} points from input file", points.len());
    points
}

fn sort_points(points: &[Point]) -> (Vec<usize>, Vec<usize>) {
    let mut sort_x: Vec<usize> = (0 .. points.len()).collect();
    sort_x.sort_by(|&i, &j| points[i].x.partial_cmp(&points[j].x).unwrap());

    let mut sort_y: Vec<usize> = (0 .. points.len()).collect();
    sort_y.sort_by(|&i, &j| points[i].y.partial_cmp(&points[j].y).unwrap());
    (sort_x, sort_y)
}

struct DiscrepancyRes {
    max_discrepancy: f64,
    rec_min_x: f64,
    rec_min_y: f64,
    rec_max_x: f64,
    rec_max_y: f64,
}

fn inner_loop_y(n: i32, low_x: f64, high_x: f64, points: &[Point],
                sub_set_y: &[usize], res: &mut DiscrepancyRes) {
    let x_dim = high_x - low_x;

    for i in 0 .. n {
        let low_y = points[sub_set_y[i as usize]].y;
        for j in i .. n {
            let high_y = points[sub_set_y[j as usize]].y;
            let y_dim = high_y - low_y;
            let count = j + 1 - i;
            let discrepancy = (points.len() as f64 * x_dim * y_dim - count as f64).abs();

            if discrepancy > res.max_discrepancy {
                *res = DiscrepancyRes {
                    max_discrepancy: discrepancy,
                    rec_min_x: low_x,
                    rec_min_y: low_y,
                    rec_max_x: high_x,
                    rec_max_y: high_y,
                };
            }
        }
    }
}

fn outer_loop_x(points: &mut [Point], sort_x: &[usize], sort_y: &[usize]) -> DiscrepancyRes {
    let mut res = DiscrepancyRes {
        max_discrepancy: 0.0,
        rec_min_x: 0.0, rec_min_y: 0.0,
        rec_max_x: 0.0, rec_max_y: 0.0,
    };

    let mut sub_set_y = vec![0; points.len()];

    for i in 0 .. points.len() {
        println!("Low X = {}", i);
        let low_x = points[sort_x[i]].x;
        for j in i .. points.len() {
            points[sort_x[j]].marked = true;
            let high_x = points[sort_x[j]].x;
            let mut n = 0;

            for k in 0 .. points.len() {
                if points[sort_y[k]].marked {
                    sub_set_y[n as usize] = sort_y[k];
                    n += 1;
                }
            }
            inner_loop_y(n, low_x, high_x, points, &sub_set_y, &mut res);
        }

        for &px in sort_x {
            points[px].marked = false;
        }
    }

    res
}

fn main() {
    let mut points = read_data();
    let (sort_x, sort_y) = sort_points(&points);
    let res = outer_loop_x(&mut points, &sort_x, &sort_y);

    println!("Maximum Discrepancy = {:.6}", res.max_discrepancy);
    println!("Rect {:.6} {:.6} {:.6} {:.6}",
             res.rec_min_x, res.rec_min_y, res.rec_max_x, res.rec_max_y);
}

Unfortunately this Rust code is a bit too much slower than the C++ version (with the standard input file the C++ code runs in 1.42 seconds, this first Rust version runs in about 2.05 seconds), probably because of the double indirection caused by the points array plus array bounds tests.

So I've rewritten the Rust code using raw pointers and lot of unsafety:

#![allow(trivial_casts)]

struct Point {
    x: f64,
    y: f64,
    marked: bool
}

fn read_data() -> Vec<Point> {
    use std::fs::File;
    use std::io::Read;

    let mut data = String::new();
    File::open("discrepancy.dat").unwrap().read_to_string(&mut data).unwrap();

    let points: Vec<Point> =
        data
        .lines()
        .map(|line| line.split(' ').map(|p| p.parse().unwrap()))
        .map(|mut parts| Point { x: parts.next().unwrap(),
                                 y: parts.next().unwrap(),
                                 marked: false })
        .collect();

    println!("Read {} points from input file", points.len());
    points
}

unsafe fn sort_points(points: &mut [Point]) -> (Vec<*mut Point>, Vec<*mut Point>) {
    let mut sort_x: Vec<*mut Point> = points.iter_mut().map(|r| r as *mut Point).collect();
    sort_x.sort_by(|&a, &b| (*a).x.partial_cmp(&(*b).x).unwrap());

    let mut sort_y: Vec<*mut Point> = points.iter_mut().map(|r| r as *mut Point).collect();
    sort_y.sort_by(|&a, &b| (*a).y.partial_cmp(&(*b).y).unwrap());
    (sort_x, sort_y)
}

struct DiscrepancyRes {
    max_discrepancy: f64,
    rec_min_x: f64,
    rec_min_y: f64,
    rec_max_x: f64,
    rec_max_y: f64,
}

unsafe fn inner_loop_y(n: i32, low_x: f64, high_x: f64, points: &[Point],
                       sub_set_y: &[*mut Point], res: &mut DiscrepancyRes) {
    let x_dim = high_x - low_x;

    for i in 0 .. n {
        let low_y = (*sub_set_y[i as usize]).y;
        for j in i .. n {
            let high_y = (*sub_set_y[j as usize]).y;
            let y_dim = high_y - low_y;
            let count = j + 1 - i;
            let discrepancy = (points.len() as f64 * x_dim * y_dim - count as f64).abs();

            if discrepancy > res.max_discrepancy {
                *res = DiscrepancyRes {
                    max_discrepancy: discrepancy,
                    rec_min_x: low_x,
                    rec_min_y: low_y,
                    rec_max_x: high_x,
                    rec_max_y: high_y,
                };
            }
        }
    }
}

unsafe fn outer_loop_x(points: &[Point],
                       sort_x: &[*mut Point], sort_y: &[*mut Point]) -> DiscrepancyRes {
    let mut res = DiscrepancyRes {
        max_discrepancy: 0.0,
        rec_min_x: 0.0, rec_min_y: 0.0,
        rec_max_x: 0.0, rec_max_y: 0.0,
    };

    let mut sub_set_y = vec![std::ptr::null::<Point>() as *mut Point; points.len()];

    for i in 0 .. points.len() {
        println!("Low X = {}", i);
        let low_x = (*sort_x[i]).x;
        for j in i .. points.len() {
            (*sort_x[j]).marked = true;
            let high_x = (*sort_x[j]).x;
            let mut n = 0;

            for k in 0 .. points.len() {
                if (*sort_y[k]).marked {
                    sub_set_y[n as usize] = sort_y[k];
                    n += 1;
                }
            }
            inner_loop_y(n, low_x, high_x, points, &sub_set_y, &mut res);
        }

        for &px in sort_x {
            (*px).marked = false;
        }
    }

    res
}

fn main() {
    let mut points = read_data();
    let res = unsafe {
        let (sort_x, sort_y) = sort_points(&mut points);
        outer_loop_x(&points, &sort_x, &sort_y)
    };

    println!("Maximum Discrepancy = {:.6}", res.max_discrepancy);
    println!("Rect {:.6} {:.6} {:.6} {:.6}",
             res.rec_min_x, res.rec_min_y, res.rec_max_x, res.rec_max_y);
}

Regarding this second Rust version, do you know if there are better ways to write these two lines of code (the #![allow(trivial_casts)] is needed by the first line)?

let mut sort_x: Vec<*mut Point> = points.iter_mut().map(|r| r as *mut Point).collect();

let mut sub_set_y = vec![std::ptr::null::<Point>() as *mut Point; points.len()];

This second Rust versions runs in about 1.85 seconds, it's slower than the C++ code, but the difference is more acceptable. How to improve this second Rust version and remove all or most unsafety? Most of the slowdown is in the inner_loop_y() function. I've asked on #IRC, and I've modified the code according to an answer by vfs:

#[derive(Clone, Copy)]
struct Point {
    x: f64,
    y: f64,
    marked: bool
}

fn read_data() -> Vec<Point> {
    use std::fs::File;
    use std::io::Read;

    let mut data = String::new();
    File::open("discrepancy.dat").unwrap().read_to_string(&mut data).unwrap();

    let points: Vec<Point> =
        data
        .lines()
        .map(|line| line.split(' ').map(|p| p.parse().unwrap()))
        .map(|mut parts| Point { x: parts.next().unwrap(),
                                 y: parts.next().unwrap(),
                                 marked: false })
        .collect();

    println!("Read {} points from input file", points.len());
    points
}

fn sort_points(mut points: Vec<Point>) -> (Vec<Point>, Vec<usize>) {
    points.sort_by(|a, b| a.x.partial_cmp(&b.x).unwrap());

    let mut sort_y: Vec<usize> = (0 .. points.len()).collect();
    sort_y.sort_by(|&i, &j| points[i].y.partial_cmp(&points[j].y).unwrap());
    (points, sort_y)
}

struct DiscrepancyRes {
    max_discrepancy: f64,
    rec_min_x: f64,
    rec_min_y: f64,
    rec_max_x: f64,
    rec_max_y: f64,
}

unsafe fn inner_loop_y(n: i32, low_x: f64, high_x: f64, points: &[Point],
                sub_set_y: &[*const Point], res: &mut DiscrepancyRes) {
    let x_dim = high_x - low_x;

    for i in 0 .. n {
        let low_y = (*sub_set_y[i as usize]).y;
        for j in i .. n {
            let high_y = (*sub_set_y[j as usize]).y;
            let y_dim = high_y - low_y;
            let count = j + 1 - i;
            let discrepancy = (points.len() as f64 * x_dim * y_dim - count as f64).abs();

            if discrepancy > res.max_discrepancy {
                *res = DiscrepancyRes {
                    max_discrepancy: discrepancy,
                    rec_min_x: low_x,
                    rec_min_y: low_y,
                    rec_max_x: high_x,
                    rec_max_y: high_y,
                };
            }
        }
    }
}

unsafe fn outer_loop_x(sort_x: &mut [Point], sort_y: &[usize]) -> DiscrepancyRes {
    let mut res = DiscrepancyRes {
        max_discrepancy: 0.0,
        rec_min_x: 0.0, rec_min_y: 0.0,
        rec_max_x: 0.0, rec_max_y: 0.0,
    };

    let mut sub_set_y: Vec<*const Point> = Vec::with_capacity(sort_x.len());

    for i in 0 .. sort_x.len() {
        println!("Low X = {}", i);
        let low_x = sort_x[i].x;
        for j in i .. sort_x.len() {
            sort_x[j].marked = true;
            let high_x = sort_x[j].x;
            let mut n = 0;
            sub_set_y.clear();

            for k in 0 .. sort_x.len() {
                if sort_x[sort_y[k]].marked {
                    sub_set_y.push(&sort_x[sort_y[k]]);
                    n += 1;
                }
            }
            inner_loop_y(n, low_x, high_x, sort_x, &sub_set_y, &mut res);
        }

        for j in 0 .. sort_x.len() {
            sort_x[j].marked = false;
        }
    }

    res
}

fn main() {
    unsafe {
        let (mut sort_x, sort_y) = sort_points(read_data());
        let res = outer_loop_x(&mut sort_x, &sort_y);

        println!("Maximum Discrepancy = {:.6}", res.max_discrepancy);
        println!("Rect {:.6} {:.6} {:.6} {:.6}",
                 res.rec_min_x, res.rec_min_y, res.rec_max_x, res.rec_max_y);
    }
}

This third Rust version contains less unsafety. The data is owned by the sorted array of Points sort_x, and I've removed the points array. The sort_y is an array of indexes, and to keep performance inner_loop_y() receives a sub_set_y that contains pointers (but now they can be const). The performance is about the same as the second Rust version, it runs in about 1.86 seconds.

Do you have suggestions for improvements and to keep (or improve) speed and remove some (or all) unsafety?

2

For kicks and giggles, have you tried using unchecked indexing of the slices? Would be interesting to see how much bounds checking adds.

3

If I replace all array accesses in the two functions inner_loop_y() and outer_loop_x() in the first Rust version that uses only indexes, its run-time is about 1.96 seconds. And doing almost the same on the second Rust program lowers its run-time to about 1.84 seconds. So array bound tests aren't the problem. I guess the problem is double indirection.

4

Well, even your unsafe version using &[*mut Point] is substantially slower than the C++ version, and that's the moral equivalent of C++. May be worthwhile to look at the generated asm.

Are you compiling C++ with clang?

5

May be worthwhile to look at the generated asm.

Right, but it's a bit complex to do. I will try to do it later.

Are you compiling C++ with clang?

I don't have Clang, I am using GCC v.6.3.0, with -Ofast. Is someone that has Clang willing to compare it with the Rust version(s)?

6

Yes, this would be interesting to compare - it might be a gcc vs clang/LLVM issue (AFAIK, gcc tends to generate better code sometimes), at least for your very unsafe version.

7

Are you using a release build (`cargo build --release && ./target/xxx" for the timing tests? If so, I'm curious how much longer the debug build takes, since that should be doing more bounds checking than the release build.

8

Are you using a release build

Yes, of course.

If so, I'm curious how much longer the debug build takes,

The release builds of the the three Rust versions run in about 51.7, 48.8 and 49.3 seconds.

9

Don't you have a discrepency?

on ccp:

        for (size_t j = i; j < points.size(); ++j)
            sortX[j]->marked = false;

on rust:

        for &px in sort_x {
            points[px].marked = false;
        }

Shouldn't it be:

       for &px in &sort_x[i..] {
            points[px].marked = false;
        }
10

And could you try with that "safe" version:

fn inner_loop_y(low_x: f64, high_x: f64, points: &[Point],
                sub_set_y: &[usize], res: &mut DiscrepancyRes) {
    let x_dim = high_x - low_x;

    for (i, &si) in sub_set_y.iter().enumerate() {
        let low_y = points[si].y;
        for (j, &sj) in sub_set_y.iter().skip(i).enumerate() {
            let high_y = points[sj].y;
            let y_dim = high_y - low_y;
            // we enumerate AFTER the skip, so we don't need to remove `i`
            let count = (j + 1) as f64;
            let discrepancy = (points.len() as f64 * x_dim * y_dim - count).abs();

            if discrepancy > res.max_discrepancy {
                *res = DiscrepancyRes {
                    max_discrepancy: discrepancy,
                    rec_min_x: low_x,
                    rec_min_y: low_y,
                    rec_max_x: high_x,
                    rec_max_y: high_y,
                };
            }
        }
    }
}

fn outer_loop_x(points: &mut [Point], sort_x: &[usize], sort_y: &[usize]) -> DiscrepancyRes {
    let mut res = DiscrepancyRes {
        max_discrepancy: 0.0,
        rec_min_x: 0.0, rec_min_y: 0.0,
        rec_max_x: 0.0, rec_max_y: 0.0,
    };

    for (i, &xi) in sort_x.iter().enumerate() {
//         println!("Low X = {}", i);
        let low_x = points[xi].x;
        for &xj in sort_x.iter().skip(i) {
            points[xj].marked = true;
            let high_x = points[xj].x;
            
            let sub_set_y = sort_y.iter()
                .filter(|&&yk| points[yk].marked)
                .cloned()
                .collect::<Vec<_>>();
            inner_loop_y(low_x, high_x, points, &sub_set_y, &mut res);
        }

        for &px in sort_x.iter().skip(i) {
            points[px].marked = false;
        }
    }

    res
}
11

Shouldn't it be:

Right, it's a porting bug... The output is the same, and the run-time seems the same, lost in the noise.

12

And could you try with that "safe" version:

Yes, it seems to work. It's the slowest of the four Rust versions, run-time about 2.09 seconds (a bit slower than the first version).

13

Sorry about that, was faster on my machine.

14

-Ofast could have something to do with it, since it turns on -ffast-math which allows GCC to make some extra assumtions about floating point values that disregards standards compliance, and thus lets it make some extra assumptions at the expense of some precision. I would assume this isn't enabled by default with rust (though LLVM does support it). Does compiling with -O3 instead make a difference?

1 Like
15

Have you tried anything to improve the speed of reading the file in? The easiest one to try would be to use a BufReader and call read_to_string() on that. I've seen that make quite a difference on large files before.

Beyond that, it would be good to get some relative numbers on how much time the Rust vs C++ spend reading in the file before rewriting readPoints() to use a streaming approach so it's more comparable to the C++.

16

Have you tried anything to improve the speed of reading the file in?

The input file is tiny, and reading it takes nothing compared to the processing time.

17

Does compiling with -O3 instead make a difference?

Yes, if I compile the C++ code (using G++ 6.3.0) with -O3 the run-time becomes even lower, about 1.36 seconds (with -Ofast it runs in about 1.41 seconds).

18

I got it down to 1.79~1.81s on my box by getting rid of marked entirely and using iterators as much as possible. There are no unsafe blocks in this version. I also left out sort() as it's own function because the lifetime bounds on sorted_y are unpleasant, but the overall result is slightly fewer lines and a little more idiomatic.

FWIW the BufReader did take off about 50ms.

#[derive(Clone, Copy)]
struct Point {
    x: f64,
    y: f64,
}

fn read_data() -> Vec<Point> {
    use std::fs::File;
    use std::io::Read;

    let mut data = String::new();
    File::open("discrepancy.dat").unwrap().read_to_string(&mut data).unwrap();

    let points: Vec<Point> =
        data
        .lines()
        .map(|line| line.split(' ').map(|p| p.parse().unwrap()))
        .map(|mut parts| Point { x: parts.next().unwrap(),
                                 y: parts.next().unwrap(),
                                })
        .collect();

    println!("Read {} points from input file", points.len());
    points
}


struct DiscrepancyRes {
    max_discrepancy: f64,
    rec_min_x: f64,
    rec_min_y: f64,
    rec_max_x: f64,
    rec_max_y: f64,
}

fn inner_loop_y(low_x: f64, high_x: f64, points: &[Point],
                sort_y: &[&Point], res: &mut DiscrepancyRes) {
    let x_dim = high_x - low_x;

    let marked = |p: &&&Point| low_x <= p.x && p.x <= high_x;

    let subset: Vec<_> = sort_y.iter().filter(marked).cloned().collect();
    for (i, p_low_y) in subset.iter().enumerate() {
        let low_y = p_low_y.y;
        for (j, p_high_y) in subset[i..].iter().enumerate() {
            let high_y = p_high_y.y;
            let y_dim = high_y - low_y;
            let count = j + 1;
            let discrepancy = (points.len() as f64 * x_dim * y_dim - count as f64).abs();

            if discrepancy > res.max_discrepancy {
                *res = DiscrepancyRes {
                    max_discrepancy: discrepancy,
                    rec_min_x: low_x,
                    rec_min_y: low_y,
                    rec_max_x: high_x,
                    rec_max_y: high_y,
                };
            }
        }
    }
}

fn outer_loop_x(sort_x: &[Point], sort_y: &[&Point]) -> DiscrepancyRes {
    let mut res = DiscrepancyRes {
        max_discrepancy: 0.0,
        rec_min_x: 0.0, rec_min_y: 0.0,
        rec_max_x: 0.0, rec_max_y: 0.0,
    };

    for (i, p_low_x) in sort_x.iter().enumerate() {
        let low_x = p_low_x.x;
        for p_high_x in &sort_x[i..] {
            let high_x = p_high_x.x;

            inner_loop_y(low_x, high_x, sort_x, &sort_y, &mut res);
        }
    }

    res
}

fn main() {
    let mut sort_x = read_data();
    sort_x.sort_by(|a, b| a.x.partial_cmp(&b.x).unwrap());

    let mut sort_y: Vec<&Point> = sort_x.iter().collect();
    sort_y.sort_by(|a, b| a.y.partial_cmp(&b.y).unwrap());


    let res = outer_loop_x(&sort_x, &sort_y);

    println!("Maximum Discrepancy = {:.6}", res.max_discrepancy);
    println!("Rect {:.6} {:.6} {:.6} {:.6}",
             res.rec_min_x, res.rec_min_y, res.rec_max_x, res.rec_max_y);
}
2 Likes
19

The code looks nice :slight_smile: On my PC with a good CPU i7 it runs in about 1.73 seconds.

20

This is derived from your code, I've pulled an allocation out of a loop:

#[derive(Copy, Clone)]
struct Point { x: f64, y: f64 }

fn read_data() -> Vec<Point> {
    use std::fs::File;
    use std::io::Read;

    let mut data = String::new();
    File::open("discrepancy.dat").unwrap().read_to_string(&mut data).unwrap();

    let points: Vec<Point> =
        data
        .lines()
        .map(|line| line.split(' ').map(|p| p.parse().unwrap()))
        .map(|mut parts| Point { x: parts.next().unwrap(),
                                 y: parts.next().unwrap() })
        .collect();

    println!("Read {} points from input file", points.len());
    points
}

struct DiscrepancyRes {
    max_discrepancy: f64,
    rec_min_x: f64,
    rec_min_y: f64,
    rec_max_x: f64,
    rec_max_y: f64,
}

fn inner_loop_y(low_x: f64, high_x: f64, points: &[Point],
                sort_y: &[&Point], subset: &mut Vec<Point>, res: &mut DiscrepancyRes) {
    let x_dim = high_x - low_x;

    subset.clear();
    let marked = |p: &&&Point| low_x <= p.x && p.x <= high_x;
    subset.extend(sort_y.iter().filter(marked).cloned());

    for (i, p_low_y) in subset.iter().enumerate() {
        let low_y = p_low_y.y;
        for (j, p_high_y) in subset[i ..].iter().enumerate() {
            let high_y = p_high_y.y;
            let y_dim = high_y - low_y;
            let count = j + 1;
            let discrepancy = (points.len() as f64 * x_dim * y_dim - count as f64).abs();

            if discrepancy > res.max_discrepancy {
                *res = DiscrepancyRes {
                    max_discrepancy: discrepancy,
                    rec_min_x: low_x,
                    rec_min_y: low_y,
                    rec_max_x: high_x,
                    rec_max_y: high_y,
                };
            }
        }
    }
}

fn outer_loop_x(sort_x: &[Point], sort_y: &[&Point]) -> DiscrepancyRes {
    let mut res = DiscrepancyRes {
        max_discrepancy: 0.0,
        rec_min_x: 0.0, rec_min_y: 0.0,
        rec_max_x: 0.0, rec_max_y: 0.0,
    };

    let mut buffer = Vec::with_capacity(sort_y.len());

    for (i, p_low_x) in sort_x.iter().enumerate() {
        let low_x = p_low_x.x;
        for p_high_x in &sort_x[i ..] {
            let high_x = p_high_x.x;
            inner_loop_y(low_x, high_x, sort_x, &sort_y, &mut buffer, &mut res);
        }
    }

    res
}

fn main() {
    let mut sort_x = read_data();
    sort_x.sort_by(|a, b| a.x.partial_cmp(&b.x).unwrap());

    let mut sort_y: Vec<&Point> = sort_x.iter().collect();
    sort_y.sort_by(|a, b| a.y.partial_cmp(&b.y).unwrap());

    let res = outer_loop_x(&sort_x, &sort_y);

    println!("Maximum Discrepancy = {:.6}", res.max_discrepancy);
    println!("Rect {:.6} {:.6} {:.6} {:.6}",
             res.rec_min_x, res.rec_min_y, res.rec_max_x, res.rec_max_y);
}

The run-time is about 1.65 seconds.

2 Likes