Closest Pair heuristics, graph with adjacency matrix in C++17

Question

I was trying to solve a problem that was briefly mentioned at the beginning of the "The Algorithm Design Manual" by Steven Skiena (Ch 1, Problem 26).

It took me some time to build a working program from the pseudocode, and I think I've got it pretty close to the described idea. However my C++ knowledge is lacking, and I'm pretty sure there must exist much easier way to achieve the goal. There is a lot of things that I doubt about, specifically:

• I have two versions of DFS-search, which seems excessive
• Four nested loops to get the pairs, is there a way to make it more human-readable? Is the complexity of that block still O(n^2)? Will I be correct if I say, that complexity of the entire solution is also O(n^2), where n - number of input points, or it's actually worse than that?
• Are there any obvious ways to make my code more clean, concise, better structured logically? Are they some well-known C++ constructions that I'm missing?
• I'm specifically interested in help, when it's possible to save lines of code without sacrificing clarity (I know it's subjective, but if there's a way to rewrite a while loop into a for loop, such that it looks clearer and takes less space, I would like to know.

I would like someone to review my code with the full rigor, and help me to improve upon it, as if my goal would be to provide a perfect C++ solution to a given problem.

The problem goes as follow:

Solution that I come up with:

#include <iostream>
#include <vector>
#include <string>
#include <cmath>

typedef std::pair<double, double> pt_t;
typedef std::vector<pt_t> pts_t;
typedef std::vector<std::vector<int>> matrix_t;

void print_point(pt_t pt) {
    std::cout << "(" << pt.first << ", " << pt.second << ")" << '\n';
}

void print_points(std::string headline, pts_t points) {
    std::cout << headline << '\n';
    std::for_each(points.begin(), points.end(), print_point);
    std::cout << "---\n";
}

void print_matrix(std::string headline, matrix_t matrix) {
    std::cout << headline << '\n';

    for (auto& row: matrix) {
        for (auto& item : row) {
            std::cout << item << ' ';
        }
        std::cout << '\n';
    }

    std::cout << "---\n";
}

void print_endpoint_pairs(std::vector<pt_t>& pairs) {
    for (auto pair : pairs) {
        std::cout << "Pair: " << pair.first << ' ' << pair.second << '\n';
    }
    std::cout << "---\n";
}

double compute_distance(const pt_t& pt1, const pt_t& pt2) {
    return std::sqrt(
        std::pow((pt1.first - pt2.first), 2) +
        std::pow((pt1.second - pt2.second), 2)
    );
}

void dfs(matrix_t& matrix, std::vector<bool>& visited, std::vector<int>& path, int v) {
    visited[v] = 1;

    path.push_back(v);

    for (int i = 0; i < matrix.size(); i++) {
        if (matrix[v][i] == 1 && !visited[i]) {
            dfs(matrix, visited, path, i);
        }
    }
}

void dfs_ep(matrix_t& matrix, std::vector<bool>& visited, std::vector<int>& path, int v) {
    visited[v] = 1;

    int connections = 0;

    for (int i = 0; i < matrix.size(); i++) {
        if (matrix[v][i] == 1) {
            connections++;
        }
    }

    // exclude points that have max number of connections
    if (connections <= 1) {
        path.push_back(v);
    }

    for (int i = 0; i < matrix.size(); i++) {
        if (matrix[v][i] == 1 && !visited[i]) {
            dfs_ep(matrix, visited, path, i);
        }
    }
}

class PlaneVector {
public:
    pts_t points{};
    matrix_t matrix;

    PlaneVector(pts_t points) :
        points(points),
        matrix(points.size(), std::vector<int>(points.size(), 0))
    {}

    matrix_t get_vertex_endpoints() {
        matrix_t chains;
        std::vector<int> chain;
        std::vector<bool> visited(points.size(), 0);

        // print_matrix("Matrix: ", matrix);

        for (int i = 0; i < points.size(); i++) {
            if (visited[i]) {
                continue;
            }

            chain.clear();

            dfs_ep(matrix, visited, chain, i);

            chains.push_back(chain);
        }

        return chains;
    }

    pts_t get_path() {
        std::vector<bool> visited(points.size(), 0);
        std::vector<int> path;
        pts_t path_points;

        dfs(matrix, visited, path, 0);

        for (int i = 0; i < path.size(); i++) {
            pt_t pt = points[path[i]];
            path_points.push_back(pt);
        }

        path_points.push_back(path_points[0]);

        return path_points;
    }

    void add_edge(int m, int n) {
        // std::cout << "Add edge: " << m << ' ' << n << '\n';
        matrix[m][n] = 1;
        matrix[n][m] = 1;
    }
};

std::vector<pt_t> get_distinct_pairs(PlaneVector& vec) {
    std::vector<pt_t> pairs{};

    matrix_t chains = vec.get_vertex_endpoints();
    // print_matrix("Endpoints: ", chains);

    // generate pairs from vertex chains endpoints
    for (int i = 0; i < chains.size() - 1; i++) {
        for (int j = i + 1; j < chains.size(); j++) {
            for (int n = 0; n < chains[i].size(); n++) {
                for (int k = 0; k < chains[j].size(); k++) {
                    pairs.push_back(std::make_pair(chains[i][n], chains[j][k]));
                }
            }
        }
    }

    return pairs;
}

pts_t closest_pair(PlaneVector& vec) {
    std::vector<pt_t> pairs = get_distinct_pairs(vec);

    while (!pairs.empty()) {
        // print_endpoint_pairs(pairs);

        double distance = std::numeric_limits<double>::max();
        int min_i = 0;
        int min_j = 0;

        for (auto pair : pairs) {
            double curr_distance = compute_distance(
                vec.points[pair.first],
                vec.points[pair.second]
            );

            if (curr_distance < distance) {
                min_i = pair.first;
                min_j = pair.second;
                distance = curr_distance;
            }
        }

        vec.add_edge(min_i, min_j);
        pairs = get_distinct_pairs(vec);
    }

    // connect two last endpoints to form a cycle
    // matrix_t chains = vec.get_vertex_endpoints();
    // vec.add_edge(chains[0][0], chains[0][1]);

    return vec.get_path();
}

int main() {
    // PlaneVector vec{{
    //     {-2, -2},
    //     {-2, 1},
    //     {1, 0},
    //     {2, -2},
    //     {2, 1},
    //     {5, 5},
    // }};

    PlaneVector vec{{
        {0.3, 0.2},
        {0.3, 0.4},
        {0.501, 0.4},
        {0.501, 0.2},
        {0.702, 0.4},
        {0.702, 0.2}
    }};

    // vec.add_edge(3, 4);
    // vec.add_edge(1, 2);
    // vec.add_edge(0, 1);
    // vec.add_edge(5, 0);

    pts_t path = closest_pair(vec);

    print_points("Points: ", vec.points);
    print_points("Path: ", path);

    return 0;
}

Answer 1

Generalizing your graph search function

The reason you had to write two versions of the graph search algorithm is that you merged the search operation with the action you want to perform on each node. You have to separate the two.

There are various approaches you could use. One is to create an iterator class that can be used to iterate over the graph in the desired order, so that you could just write something like:

for (auto v: dfs(matrix)) {
    path.push_back(v);
}

Alternatively, you can write a function that takes a function object as a parameter, and applies it on each node that it finds in the desired order. You also want to avoid having to pass visited and v as a parameter to dfs(), since those variables are just internal details of the DFS algorithm, you should not expose that.

static void dfs_impl(const matrix_t &matrix, std::function<void(int)> &func, static void dfs_impl(const matrix_t &matrix, const std::function<void(int)> &func, std::vector<bool> &visited, int v) {
    visited[v] = true;

    func(v);

    for (int i = 0; i < matrix.size(); ++i) {
        if (matrix[v][i] && !visited[i]) {
            dfs_impl(matrix, func, visited, i);
        }
    }
}

void dfs2(const matrix_t &matrix, int root, const std::function<void(int)> &func) {
    std::vector<bool> visited(matrix.size());
    dfs_impl(matrix, func, visited, root);
}

Now you can call it like so:

pts_t get_path() const {
    pts_t path_points;

    dfs(matrix, 0, [&](int v){ path_points.push_back(points[v]); });

    path_points.push_back(path_points.front());
    return path_points;
}

And instead of calling dfs_ep(), you can write the following:

matrix_t get_vertex_endpoints() const {
    matrix_t chains;
    std::vector<bool> visited(points.size());

    for (int i = 0; i < points.size(); i++) {
        if (visited[i]) {
            continue;
        }

        std::vector<int> chain;

        dfs(matrix, i, [&](int v){
            visited[v] = true;

            if (std::count(matrix[v].begin, matrix[v].end, 1) <= 1) {
                chain.push_back(v);
            }
        });

        chains.push_back(chain);
    }

    return chains;
}

Note that here we had to keep a local vector visited. You could make it so you still pass a reference to visited to the function dfs(), but I find this is not as clean. Another approach is to have dfs() return an iterator to the next unvisited node:

int dfs(...) {
    std::vector<bool> visited(matrix.size());
    dfs_impl(matrix, func, visited, root);
    return std::find(visited.begin() + root, visited.end(), false) - visited.begin();
}

In that case, you can rewrite get_vertex_endpoints() like so:

matrix_t get_vertex_endpoints() const {
    ...
    for (int i = 0; i < points.size();) {
        ...
        i = dfs(matrix, i, [&](int v){
            ...

Nesting for-loops

It is probably possible to make the four nested for-loops in get_distinct_pairs() look better. You could make a class that allows iteration over pairs, and use some kind of Cartesian product iterator from existing libraries, and use C++17 structured bindings to make the for-loops look approximately like this:

for (auto [chain1, chain2]: pairs(chains)) {
    for (auto [vertex1, vertex2]: cartesian_product(chain1, chain2) {
        pairs.push_back({vertex1, vertex2});
    }
}

However, those functions are not in the standard library, so to be portable you'd have to implement them yourself. I don't think four nested loops is bad here, the comment explains what you are going to do.

The complexity is still just O(n^2).

Other ways to make the code more readable

There are lots of functions in the standard library that can help you. I already shown a few example above, where I used std::count() and std::find() to remove manual loops. Not only does it make the code shorter, it also expresses intent explicitly.

There's also some places where you can use auto, structured bindings and so on to reduce the amount of code without hurting readability. I'll mention some more specific things that can be improved below.

Use std::hypot()

To compute the distance between two 2D points, you can make use of std::hypot():

double compute_distance(const pt_t& pt1, const pt_t& pt2) {
    return std::hypot(pt1.first - pt2.first, pt1.second - pt2.second);
}

Write std::ostream formatters instead of print() functions

Instead of writing print_point(pt), wouldn't it be nicer to be able to write std::cout << pt << '\n'? You can do this by converting your printing functions to overload the <<-operator of std::ostream, like so:

std::ostream &operator<<(std::ostream &o, const pt_t &pt) {
    return o << "(" << pt.first << ", " << pt.second << ")";
}

Apart from printing your own objects in a more idiomatic way, it's now also much more generic, and allows you to print to files, stringstreams, and everything else that is a std::ostream.

Use const where appropariate

Any time a function takes a pointer or reference parameter, and does not modify it, you should mark it as const, so the compiler can better optimize your code, and can give an error if you accidentily do modify it.

Also, class member functions that do not modify any of the member variables should also be marked const.