Distance between two different sets of points

I see a number of things that may help you improve your program. I'll start with the more superficial and progress to more substantive suggestions.

Don't abuse using namespace std

Putting using namespace std at the top of every program is a bad habit that you'd do well to avoid. Know when to use it and when not to (as when writing include headers). In this particular case, I happen to think it's not terrible but it also doesn't help much.

Make sure you have all required #includes

The code uses sqrt but doesn't #include <cmath>. Also, remember that sqrt is in the std namespace, so the fully qualified name is std::sqrt. Additionally, rand is used but no <cstdlib>

Be careful with signed and unsigned

In the current code, the loop integers i and j and k are signed int values, but they're being compared with unsigned quantities x1.size() and x2.size(). Better would be to declare them all as unsigned or size_t.

Use const where practical

The passed vectors should not be altered (and are not altered) by the processing functions. Indicate that fact (and maybe even get a small speed boost) by specifying const std::vector& as the function arguments.

Use "range for" and simplify your code

If you're using a C++11 compliant compiler, the use of "range for" can simplify your code. For example, the code currently contains this:

for (unsigned i = 0; i < res.size(); ++i) if (res[i] < cutoff) res2.push_back(res[i]); 

It could be this:

for (auto dist: res) if (dist < cutoff) res2.push_back(dist); 

Better would be to eliminate it completely, but we'll get to that later.

Don't use std::endl if you don't really need it

The difference between std::endl and '\n' is that '\n' just emits a newline character, while std::endl actually flushes the stream. This can be time-consuming in a program with a lot of I/O and is rarely actually needed. It's best to only use std::endl when you have some good reason to flush the stream and it's not very often needed for simple programs such as this one. Avoiding the habit of using std::endl when '\n' will do will pay dividends in the future as you write more complex programs with more I/O and where performance needs to be maximized.

Don't store data you don't need

Right now, the code attempts to calculate the comprehensive list of distances and puts them into a large vector. On my machine, a double is 8 bytes. With 50000 and 25000 points in the two sets, this means the resulting distance vector is 50000 * 25000 * 8 = 10,000,000,000 bytes or 9.3 GiB. That's a huge amount of data, which the program then mostly discards and then reduces to a single number. (On my machine, that number is always -nan which suggests another problem, but more on that later.)

Use classes more effectively

The description talks about points, but the program actually uses vectors of coordinates instead. I'd recommend creating a templated Point2D object like this:

template <typename T> class Point2D { public: Point2D(T x, T y) : x_(x), y_(y) {} T dist(const Point2D& other) const { return std::sqrt(sqdist(other)); } T sqdist(const Point2D& other) const { const auto dx{x_ - other.x_}; const auto dy{y_ - other.y_}; return dx*dx + dy*dy; } private: T x_, y_; }; 

Consider using a better random number generator

Right now, the program is using the old C-style rand which is not a very good random number generator. The intent is apparently to create numbers in the range [0, RAND_MAX] where RAND_MAX is implementation defined. On my machine, RAND_MAX = 2147483647 but on yours I suspect it must be 32768. Otherwise you would also be getting -nan for the result or you'd have chosen a different cutoff value. So rather than relying on implementation defined values and a poor random number generator, if you are using a compiler that supports at least C++11, consider using a better random number generator. In particular, instead of rand, you might want to look at std::uniform_real_distribution and friends in the <random> header. Here's one way to do it:

double newrand() { static std::random_device rd; static std::mt19937 gen(rd()); static std::uniform_real_distribution<double> dis(0, 32768); return dis(gen); } 

Avoid computationally costly operations

The distance between every pair of points is calculated which means millions of calls to std::sqrt, but this is not really needed since it would be sufficient to calculate the squared distance instead and then compare with the squared cutoff value. We can go a bit further and test the square of each coordinate pair. That is, if \$(\Delta x)^2 \ge t^2\$, or \$(\Delta y)^2 \ge t^2\$, then there's no point in doing further calculations.

Results

When I first ran the program, it takes 30.9 seconds just for the all_dist call and a total of 55.1 seconds to finally print the answer. After applying all of the suggestions above the program runs in 0.95 seconds. Here's the revised code:

#include <vector> #include <algorithm> #include <iostream> #include <ctime> #include <cmath> #include <random> double newrand() { static std::random_device rd; static std::mt19937 gen(rd()); static std::uniform_real_distribution<double> dis(0, 32768); return dis(gen); } template <typename T> class Point2D { public: Point2D() : x_(newrand()), y_(newrand()) {} Point2D(T x, T y) : x_(x), y_(y) {} T dist(const Point2D& other) const { return std::sqrt(sqdist(other)); } T sqdist(const Point2D& other) const { const auto dx{x_ - other.x_}; const auto dy{y_ - other.y_}; return dx*dx + dy*dy; } T sqdist_thr(const Point2D& other, T threshold) const { const auto dx{x_ - other.x_}; const auto dy{y_ - other.y_}; if (dx < threshold && dy < threshold) return dx*dx + dy*dy; return threshold; } private: T x_, y_; }; double avgdist(const std::vector<Point2D<double>>& a, const std::vector<Point2D<double>>& b, double threshold) { const double threshold2 = threshold*threshold; double sum{0}; unsigned count{0}; for (const auto &one: a) { for (const auto &two: b) { auto d2 = one.sqdist_thr(two, threshold2); if (d2 < threshold2) { ++count; sum += std::sqrt(d2); } } } return sum/count; } int main() { constexpr int num_el = 50000; constexpr double cutoff = 200; const std::vector<Point2D<double>> a(num_el); const std::vector<Point2D<double>> b(num_el/2); clock_t begin_time = clock(); std::cout << avgdist(a, b, cutoff) << '\n'; std::cout << float(clock() - begin_time) / CLOCKS_PER_SEC << '\n'; }