Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Math/Random and its tests implementation #432

Open
wants to merge 11 commits into
base: master
Choose a base branch
from
Prev Previous commit
Next Next commit
public functions
sariug committed Apr 11, 2020
commit c5e5804d730fe42621edbe39c418348acfbce80f
67 changes: 32 additions & 35 deletions src/Magnum/Math/Random.h
Original file line number Diff line number Diff line change
@@ -17,78 +17,75 @@ namespace Magnum
namespace Math
{

namespace Implementation
namespace Random
{
static std::seed_seq seeds{{
static_cast<std::uintmax_t>(std::random_device{}()),
static_cast<std::uintmax_t>(std::chrono::steady_clock::now()
.time_since_epoch()
.count()),
}};

class RandomGenerator
{
public:
RandomGenerator() = delete;

RandomGenerator()
{
std::seed_seq seeds{{
static_cast<std::uintmax_t>(std::random_device{}()),
static_cast<std::uintmax_t>(std::chrono::steady_clock::now()
.time_since_epoch()
.count()),
}};
g = std::mt19937{seeds};
};
template <typename T>
static typename std::enable_if<std::is_same<Int, T>::value, T>::type
typename std::enable_if<std::is_same<Int, T>::value, T>::type
generate(T start = -Magnum::Math::Constants<T>::inf(),
T end = Magnum::Math::Constants<T>::inf())
{
return std::uniform_int_distribution<T>{start, end}(generator());
return std::uniform_int_distribution<T>{start, end}(g);
}

template <typename T>
static typename std::enable_if<std::is_same<Float, T>::value, T>::type
typename std::enable_if<std::is_same<Float, T>::value, T>::type
generate(T start = -Magnum::Math::Constants<T>::inf(),
T end = Magnum::Math::Constants<T>::inf())
{
return std::uniform_real_distribution<T>{start, end}(generator());
return std::uniform_real_distribution<T>{start, end}(g);
}

public:
static std::mt19937 &generator()
{
static std::mt19937 g{seeds};
return g;
}
private:
// namespace Implementation
std::mt19937 g;
};
} // namespace Implementation
namespace Random
{

template <class T = Float>
T randomScalar(T begin = 0.0f, T end = 1.0f)
T randomScalar(RandomGenerator &g, T begin = 0.0f, T end = 1.0f)
{
return Implementation::RandomGenerator::generate(static_cast<T>(begin),
static_cast<T>(end));

return g.generate(static_cast<T>(begin),
static_cast<T>(end));
}

template <class T = Float>
Vector2<T> randomUnitVector2()
Vector2<T> randomUnitVector2(RandomGenerator &g)
{
auto a = Implementation::RandomGenerator::generate(0.0f, 2 * Math::Constants<T>::pi());
auto a = g.generate(0.0f, 2 * Math::Constants<T>::pi());
return {std::cos(a), std::sin(a)};
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I guess this is where things become hard 😅

I actually have no idea here -- will the distribution be still uniform if sin/cos is used? I guess it will? Ideally this would be without the extra overhead of trig functions, but I don't have any idea if that's doable ... in this thread on SO they use a Gaussian distribution as an input, but .. ¯\_(ツ)_/¯

If you have better references than me, mention them here please :)

Copy link
Author

@sariug sariug Apr 8, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For this case sin and cos shall not be get affected ? (sin (+ + - - )[50%], cos (+ - - + )[50%] )
So this shall be as accurate as the generator itself ?

Can you also comment about what I read here ? It looks accurate according to this. The main discussion seems like "getting 3 random numbers and normalizing" vs "2 numbers(theta and height) and calculating". The latter seems more accurate, which is sin/cos.

Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes that's a great reference, thanks -- a link to it should go in the documentation. I didn't absorb it fully yet, but yeah it sounds like a good proof.

getting 3 random numbers and normalizing

this is what you do for quaternions tho .. can the two-/three-dimensional case be extended for those as well?

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Quaternion is link here.

Added both to the code :)

}

template <class T = Float>
Vector3<T> randomUnitVector3()
Vector3<T> randomUnitVector3(RandomGenerator &g)
{
// Better to have it "theta" and "z" than three random numbers.
// https://mathworld.wolfram.com/SpherePointPicking.html
auto a = Implementation::RandomGenerator::generate(0.0f, 2 * Math::Constants<T>::pi());
auto z = randomScalar(-1.0f, -1.0f);
auto a = g.generate(0.0f, 2 * Math::Constants<T>::pi());
auto z = randomScalar(g, -1.0f, -1.0f);
auto r = sqrt<T>(1 - z * z);
return {r * std::cos(a), r * std::sin(a), z};
}

template <class T = Float>
Quaternion<T> randomRotation()
Quaternion<T> randomRotation(RandomGenerator &g)
{
//http://planning.cs.uiuc.edu/node198.html
auto u = randomScalar();
auto v = 2 * Math::Constants<T>::pi() * randomScalar();
auto w = 2 * Math::Constants<T>::pi() * randomScalar();
auto u = randomScalar(g);
auto v = 2 * Math::Constants<T>::pi() * randomScalar(g);
auto w = 2 * Math::Constants<T>::pi() * randomScalar(g);
return Quaternion<T>({sqrt<T>(1 - u) * std::sin(v),
sqrt<T>(1 - u) * std::cos(v),
sqrt<T>(u) * std::sin(w)},
22 changes: 15 additions & 7 deletions src/Magnum/Math/Test/RandomTest.cpp
Original file line number Diff line number Diff line change
@@ -23,6 +23,7 @@ struct RandomTest : Corrade::TestSuite::Tester
void unitVector3();
void randomRotation();
void randomDiceChiSquare();

};

typedef Vector<2, Float> Vector2;
@@ -43,29 +44,36 @@ RandomTest::RandomTest()

void RandomTest::randScalar()
{
CORRADE_COMPARE_AS(Math::Random::randomScalar<Float>(-1.0, 1.0), 1.0f, Corrade::TestSuite::Compare::LessOrEqual);
CORRADE_COMPARE_AS(Math::Random::randomScalar<Float>(-1.0, 1.0), -1.0f, Corrade::TestSuite::Compare::GreaterOrEqual);
Math::Random::RandomGenerator g;
CORRADE_COMPARE_AS(Math::Random::randomScalar<Float>(g, -1.0, 1.0), 1.0f, Corrade::TestSuite::Compare::LessOrEqual);
CORRADE_COMPARE_AS(Math::Random::randomScalar<Float>(g, -1.0, 1.0), -1.0f, Corrade::TestSuite::Compare::GreaterOrEqual);
}

void RandomTest::unitVector2()
{
CORRADE_COMPARE((Math::Random::randomUnitVector2()).length(), 1.0f);
Math::Random::RandomGenerator g;
CORRADE_COMPARE((Math::Random::randomUnitVector2(g)).length(), 1.0f);
}
void RandomTest::unitVector3()
{
CORRADE_COMPARE((Math::Random::randomUnitVector3()).length(), 1.0f);
Math::Random::RandomGenerator g;

CORRADE_COMPARE((Math::Random::randomUnitVector3(g)).length(), 1.0f);
}

void RandomTest::randomRotation()
{
CORRADE_COMPARE(Math::Random::randomRotation().length(), 1.0f);
Math::Random::RandomGenerator g;

CORRADE_COMPARE(Math::Random::randomRotation(g).length(), 1.0f);
}

void RandomTest::randomDiceChiSquare()
{
// A step by step explanation
// https://rpg.stackexchange.com/questions/70802/how-can-i-test-whether-a-die-is-fair

Math::Random::RandomGenerator g;

int error_count = 0; // We have 1 chance to over shoot. Thats why no repeated test.

const Int dice_side = 20;
@@ -77,7 +85,7 @@ void RandomTest::randomDiceChiSquare()

std::vector<Int> faces(dice_side, 0);
for (std::size_t i = 0; i < expected * dice_side; i++)
faces[Math::Random::randomScalar<Int>(0, dice_side - 1)]++;
faces[Math::Random::randomScalar<Int>(g, 0, dice_side - 1)]++;
std::vector<Int> residual(dice_side, 0);
for (std::size_t i = 0; i < dice_side; i++)
residual[i] = Float(pow((faces[i] - expected), 2)) / expected;