-1

I have a function to generate a (pseudo) random walk on a square lattice where the walk should not breach the boundaries of this square, full function below:

/**
* @brief Performs a single random walk returning the final distance from the origin
*
* Completes a random walk on a square lattice using the mersenne twister engine based pseudo-random
* number-generator (PRNG). The walk will not breach the boundaries of the square size provided to
* the function. The random walk starts at the origin and ends after some parameterised number of steps.
* Position co-ordinates of the walk for each iteration are sent to an output file.
*
* @param squareSideLength Length of square lattice side
* @param steps Number of steps to compute random walk up to
* @param engine Mersenne Twister engine typedef (used for generating random numbers locally)
* @param distribution Default distribution of random walk
* @param outputFile [Default nullptr] Pointer to file to write co-ordinate data of random walk to
* @return final distance of the particle from the origin
*/
double randomWalkSquareLattice(int squareSideLength, int steps, std::mt19937& engine, std::uniform_real_distribution<double>& distribution, std::ofstream* outputFile = nullptr) {

    // store the half-length of the square lattice
    const int halfSquareLength = squareSideLength / 2;

    // initialise co-ordinates to the origin
    double positionX = 0.0;
    double positionY = 0.0;

    // assign the default distribution to distDefault
    std::uniform_real_distribution<double> distDefault = distribution;

    // loop over a number of iterations given by the steps parameter
    for (int i = 0; i < steps; i++) {

        std::cout << positionX << "\t" << positionY << std::endl;

        // if the x-position of the particle is >= to positive
        // half square lattice length then generate decremental 
        // random number (avoiding breaching the boundary)
        if (positionX >= halfSquareLength) {
            double offset = positionX - halfSquareLength;
            std::cout << std::endl << offset << std::endl;
            std::uniform_real_distribution<double> distOffset(-offset, -1.0);
            positionX += distOffset(engine);
        }

        // else if the x-position of the particle is <= to negative
        // half square lattice length then generate incremental random
        // number (avoiding breaching the boundary)
        else if (positionX <= -halfSquareLength) {
            double offset = std::abs(positionX + halfSquareLength);
            std::cout << std::endl << offset << std::endl;
            std::uniform_real_distribution<double> distOffset(offset, 1.0);
            positionX += distOffset(engine);
        }

        // else (in case where x-position of particle is not touching 
        // the lattice boundary) generate default random number
        else {
            positionX += distDefault(engine);
        }

        // if the y-position of the particle is >= to positive
        // half square lattice length then generate decremental 
        // random number (avoiding breaching the boundary)
        if (positionY >= halfSquareLength) {
            double offset = positionY - halfSquareLength;
            std::cout << std::endl << offset << std::endl;
            std::uniform_real_distribution<double> distOffset(-offset, -1.0);
            positionY += distOffset(engine);
        }

        // else if the y-position of the particle is <= to negative
        // half square lattice length then generate incremental 
        // random number (avoiding breaching the boundary)
        else if (positionY <= -halfSquareLength) {
            double offset = std::abs(positionY + halfSquareLength);
            std::cout << std::endl << offset << std::endl;
            std::uniform_real_distribution<double> distOffset(offset, 1.0);
            positionY += distOffset(engine);
        }

        // else (in case where y-position of particle is not touching
        // the lattice boundary) generate default random number
        else {
            positionY += distDefault(engine);
        }

        // if an outputFile is supplied to the function, then write data to it
        if (outputFile != nullptr) {
            *outputFile << positionX << "\t" << positionY << std::endl;
        }

    }

    // compute final distance of particle from origin
    double endDistance = std::sqrt(positionX*positionX + positionY*positionY);

    return endDistance;

}

Where the conditionals seen in the method prevent the walk exiting the boundaries. However, when this is called with a sufficient number of steps (so that any one of these conditionals is executed) I get an error saying:

invalid min and max arguments for uniform_real

Note that the dist I send to this function is:

std::uniform_real_distribution<double> dist(-1.0,1.0);

And so (as you can see from the values printed to the terminal) the issue is not that the offset will ever be larger than the max value given to the distOffset in any of the conditional cases.

Is the issue that I cannot give u_r_d a double value of arbitrary precision? Or is something else at play here that I am missing?

Edit: I should add that these are the values used in main():

int main(void) {
    std::uniform_real_distribution<double> dist(-1.0, 1.0);

    std::random_device randDevice;

    std::mt19937 engine(randDevice());
    //std::cout << dist(engine) << std::endl;
    // Dimensions of Square Lattice
    const int squareLength = 100;

    // Number of Steps in Random Walk
    const int nSteps = 10000;

    randomWalkSquareLattice(squareLength, nSteps, engine, dist);
}
sjrowlinson
  • 3,297
  • 1
  • 18
  • 35
  • 1
    Time to use a debugger. There are lots of places where you construct a `uniform_real_distribution` with a variable upper or lower bound. – MicroVirus Nov 23 '15 at 18:43
  • Why do you create a new `uniform_real_distribution` every time? You could simply create one `uniform_real_distribution` in `(0, 1)`, obtain a sample, and multiply it by the desired offset. The result would be a uniform integer in the range `(0, desired_offset)`. – Escualo Nov 23 '15 at 18:48
  • Off topic: `std::mt19937 engine(randDevice());` beware of this. You've done it right, but... http://stackoverflow.com/questions/18880654/why-do-i-get-same-sequence-for-everyrun-with-stdrandom-device-with-mingw-gcc4 – user4581301 Nov 23 '15 at 18:51
  • @Escualo I don't think that will work as I need to generate a random double with a minimum value of the offset (i.e. the amount by which the quantity has breached a given boundary) and a maximum value of (e.g.) 1.0. Using the method you outlined would, in general, allow the random walk to travel outside the boundary or get stuck on it. I suppose, however, I could alter your method by adding the offset to the [0,1] interval and so I'll try that instead. – sjrowlinson Nov 23 '15 at 19:25
  • What I'm trying to say is that having a random number in [0, 1] lets you generate a number in an arbitrary interval: `y = y_lb + x * (y_ub - y_lb)`, where `y` is your new random number, `y_lb` is the desired lower bound, `y_ub` is the desired upper bound, and `x` is a random number in `[0, 1]`. How you pick `y_lb` and `y_ub` is up to you. – Escualo Nov 23 '15 at 19:31
  • I have managed to achieve the desired result from using the method I mentioned above, thanks for the assistance. I can post the code in an edit to the OP if either of you would like me to. Also, I assume that using engine(randDevice()) will not be a problem? It seems to generate different random data each run for me anyway. – sjrowlinson Nov 23 '15 at 19:46
  • Different tool-chains, different results. I'd leave a note in the code to help you out if you port to MinGw and they haven't sorted out the bug, though. Could save future time. Also, don't edit the question. Make an answer. – user4581301 Nov 23 '15 at 19:48

2 Answers2

0

uniform_real_distribution(a,b); requires that a ≤ b.

If positionX == halfSquareLength, then,

double offset = positionX - halfSquareLength;

is the same as saying

double offset = positionX - positionX;

and offset will be zero.

This results in

std::uniform_real_distribution<double> distOffset(-0.0, -1.0);

and violates a ≤ b.

user4581301
  • 33,082
  • 7
  • 33
  • 54
0

Here is the solution I came up with, seems to work for all test cases so far:

/**
* @brief Performs a single random walk returning the final distance from the origin
*
* Completes a random walk on a square lattice using the mersenne twister engine based pseudo-random
* number-generator (PRNG). The walk will not breach the boundaries of the square size provided to
* the function. The random walk starts at the origin and ends after some parameterised number of steps.
* Position co-ordinates of the walk for each iteration are sent to an output file.
*
* @param squareSideLength Length of square lattice side
* @param steps Number of steps to compute random walk up to
* @param engine Mersenne Twister engine typedef (used for generating random numbers locally)
* @param distribution Default distribution of random walk
* @param outputFile [Default nullptr] Pointer to file to write co-ordinate data of random walk to
* @return final distance of the particle from the origin
*/
double randomWalkSquareLattice(int squareSideLength, int steps, std::mt19937& engine, std::uniform_real_distribution<double>& distribution, std::ofstream* outputFile = nullptr) {

// store the half-length of the square lattice
const int halfSquareLength = squareSideLength / 2;

// initialise co-ordinates to the origin
double positionX = 0.0;
double positionY = 0.0;

// assign the default distribution to distDefault
std::uniform_real_distribution<double> distDefault = distribution;

std::uniform_real_distribution<double> distBound(0.0, 1.0);

double oS;

// loop over a number of iterations given by the steps parameter
for (int i = 0; i < steps; i++) {

    //std::cout << positionX << "\t" << positionY << std::endl;

    positionX += distDefault(engine);
    positionY += distDefault(engine);

    // if the x-position of the particle is >= to positive
    // half square lattice length then generate decremental 
    // random number (avoiding breaching the boundary)
    if (positionX >= halfSquareLength) {
        oS = distBound(engine);
        double offset = positionX - halfSquareLength;
        double desiredOffset = -(oS + offset);

        if (desiredOffset < -1.0) {
            double offsetFromNegUnity = desiredOffset + 1.0;
            desiredOffset -= offsetFromNegUnity;
        }

        positionX += desiredOffset;
    }

    // else if the x-position of the particle is <= to negative
    // half square lattice length then generate incremental random
    // number (avoiding breaching the boundary)
    else if (positionX <= -halfSquareLength) {
        oS = distBound(engine);
        double offset = std::abs(positionX + halfSquareLength);
        double desiredOffset = offset+oS;

        if (desiredOffset > 1.0) {
            double offsetFromUnity = desiredOffset - 1.0;
            desiredOffset -= offsetFromUnity;
        }

        positionX += desiredOffset;
    }

    // if the y-position of the particle is >= to positive
    // half square lattice length then generate decremental 
    // random number (avoiding breaching the boundary)
    if (positionY >= halfSquareLength) {
        oS = distBound(engine);
        double offset = positionY - halfSquareLength;
        double desiredOffset = -(offset+oS);

        if (desiredOffset < -1.0) {
            double offsetFromNegUnity = desiredOffset + 1.0;
            desiredOffset -= offsetFromNegUnity;
        }

        positionY += desiredOffset;
    }

    // else if the y-position of the particle is <= to negative
    // half square lattice length then generate incremental 
    // random number (avoiding breaching the boundary)
    else if (positionY <= -halfSquareLength) {
        oS = distBound(engine);
        double offset = std::abs(positionY + halfSquareLength);
        double desiredOffset = offset+oS;

        if (desiredOffset > 1.0) {
            double offsetFromUnity = desiredOffset - 1.0;
            desiredOffset -= offsetFromUnity;
        }

        positionY += desiredOffset;
    }

    // if an outputFile is supplied to the function, then write data to it
    if (outputFile != nullptr) {
        *outputFile << positionX << "\t" << positionY << std::endl;
    }

}

// compute final distance of particle from origin
double endDistance = std::sqrt(positionX*positionX + positionY*positionY);

return endDistance;

}

Here, an offset was generated randomly on the interval (0,1) and the difference from the boundary by which a x or y position breached was added to this offset to create a double value which would have a minimum of this breaching difference and (after an additional nested conditional check) a maximum of 1.0 (or -1.0 for opposite boundary).

sjrowlinson
  • 3,297
  • 1
  • 18
  • 35