-1

I want to create an MxN array (M particles in N dimensional space) filled with random numbers within an upper and lower boundary. I have a working python code that looks something like this:

# upper_bound/lower_bound are arrays of shape (dim,)
positions = np.random.rand(num_particle,dim)*(upper_bound-lower_bound)+lower_bound

Each row represents a particle, and each column represents a dimension in the problem space. So the upper_bound and lower_bound applies to each column. Now I want to translate the above code to c++, and I have something like this:

#include <iostream>
#include <vector>
#include <random>
#include <algorithm>
#include <ctime>

typedef std::vector<double> vect1d;

std::vector<vect1d> positions;

for (int i=0; i<num_particle; i++){
    std::mt19937_64 generator(static_cast<std::mt19937::result_type>(time(0)));
    std::uniform_real_distribution<double> distribution(0,1);
    vect1d pos(dimension);
    std::generate(pos.begin(),pos.end(),distribution(generator));
    positions[i] = pos;
    }

My problems:

  1. It gives error regarding the generator, so I'm not sure if I set it properly. I'm also not sure how to use the std::generator. I'm trying it as I've looked at other similar posts and it seems that it allows me to generate more than one random number at a time, so I don't have to run it MxN times for each element. Is this true and how to correctly use it?

  2. In python I can just vectorization and broadcasting to manipulate the numpy array. What's the most 'vectorized' way to do it in c++?

  3. The above (incorrect) code only creates random number between 0 and 1, but how to incorporate the lower_bound and upper_bound as in the python version? I understand that I can change the values inside distribution(0,1), but the problem is the limits can be different for each dimension (so each column can have different valid range), so what's the most efficient way to generate random number, taking into account the range for each dimension?

Thanks

Physicist
  • 2,848
  • 8
  • 33
  • 62
  • 1
    Do you set size of `positions`? If it is empty, assigning to `positions[i]` is UB. – Evg Sep 01 '18 at 13:28
  • Are you on Visual Studio? – sandthorn Sep 01 '18 at 13:31
  • 1
    How to use a generator: `std::generate(pos.begin(), pos.end(), [&]() { return distribution(generator); });` – Evg Sep 01 '18 at 13:35
  • For future reference it's best to split these sorts of multi-questions into separate ones. Increase the chance you're getting some answers. – Horia Coman Sep 01 '18 at 13:39
  • This may be helpful with `3.` https://stackoverflow.com/questions/52021277/choose-random-number-distribution-at-compile-time/52021446#52021446 Need to use floating point arguments to get floating point results aka `random_number(0.0, 1.0)`. – Galik Sep 01 '18 at 13:52
  • Using a library like Eigen would simplify a few things. – Marc Glisse Sep 04 '18 at 07:58

2 Answers2

1

I'll address them in random order:

  • 3.You have several options - using one generator per row, created like distribution(row_lower_limit, row_upper_limit). Should be cheap enough to not cause issues. If you want to reuse the same generator, just do something like row_lower_limit + distribution(generator) * (row_upper_limit - row_lower_limit). The distribution is in both cases U[row_lower_limit, row_upper_limit].
  • 2.The vectorization came from the numpy library, not from Python itself. It provided some nice UX at most. C++ doesn't have an equivalent library to numpy (though there's a lot of libraries for it as well - just nothing so univeral). You wouldn't be wrong by doing two nested fors. You'd perhaps be better served by just declaring a NxM array rather than a vector, like here.
  • 1.Not sure how to help with the problem since we don't know the error. The cplusplus.com reference has an example of how to initialize this with reference to a random_device.
Horia Coman
  • 8,681
  • 2
  • 23
  • 25
1

First of all, you're doing more work than you need to with your Python version, just use:

np.random.uniform(lower_bound, upper_bound, size=(num_particle, dim))

In your C++ attempt, the line

std::generate(pos.begin(),pos.end(),distribution(generator));

Is incorrect as the third argument must be a function not a value. A reasonable C++ equivalent would be:

using RandomVector = std::vector<double>;
using RandomMatrix = std::vector<RandomVector>;

template <typename Generator=std::mt19937_64>
RandomMatrix&
fill_uniform(const double low, const double high, RandomMatrix& result)
{
    Generator gen {static_cast<typename Generator::result_type>(time(0))};
    std::uniform_real_distribution<double> dist {low, high};
    for (auto& col : result) {
        std::generate(std::begin(col), std::end(col), [&] () { return dist(gen); });
    }
    return result;
}

template <typename Generator=std::mt19937_64>
RandomMatrix
generate_uniform(const double low, const double high,
                 const std::size_t ncols, const std::size_t nrows)
{
    RandomMatrix result(ncols, RandomVector(nrows));
    return fill_uniform<Generator>(low, high, result);
}

int main()
{
    auto m = generate_uniform(2, 11, 2, 3);
    for (const auto& col : m) {
        for (const auto& v : col) {
            std::cout << v << " ";
        }
        std::cout << '\n';
    }
}

You could generalise this to generate arbitrary dimension tensors (like the NumPy version) without too much work.

Daniel
  • 8,179
  • 6
  • 31
  • 56