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
1 change: 1 addition & 0 deletions src/Magnum/Math/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -49,6 +49,7 @@ set(MagnumMath_HEADERS
Quaternion.h
Packing.h
PackingBatch.h
Random.h
Range.h
RectangularMatrix.h
StrictWeakOrdering.h
Expand Down
100 changes: 100 additions & 0 deletions src/Magnum/Math/Random.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,100 @@
#ifndef Magnum_Math_Random_h
#define Magnum_Math_Random_h

// TO DO Licence things.

#include <random>
#include <chrono>
#include "Magnum/Types.h"
#include "Magnum/Math/Constants.h"
#include "Magnum/Math/Vector2.h"
#include "Magnum/Math/Vector3.h"
#include "Magnum/Math/Quaternion.h"
#include "Magnum/Math/Functions.h"

namespace Magnum
{
namespace Math
{

namespace Random
{
class RandomGenerator
{
public:
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>
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}(g);
}

template <typename T>
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}(g);
}

private:
// namespace Implementation
std::mt19937 g;
};

template <class T = Float>
T randomScalar(RandomGenerator &g, T begin = 0.0f, T end = 1.0f)
{

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

template <class T = Float>
Vector2<T> randomUnitVector2(RandomGenerator &g)
{
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(RandomGenerator &g)
{
// Better to have it "theta" and "z" than three random numbers.
// https://mathworld.wolfram.com/SpherePointPicking.html
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(RandomGenerator &g)
{
//http://planning.cs.uiuc.edu/node198.html
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)},
sqrt<T>(u) * std::cos(w));
}

} // namespace Random

} // namespace Math
} // namespace Magnum

#endif
6 changes: 6 additions & 0 deletions src/Magnum/Math/Test/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -70,6 +70,8 @@ corrade_add_test(MathStrictWeakOrderingTest StrictWeakOrderingTest.cpp LIBRARIES

corrade_add_test(MathMatrixBenchmark MatrixBenchmark.cpp LIBRARIES MagnumMathTestLib)

corrade_add_test(MathRandomTest RandomTest.cpp LIBRARIES MagnumMathTestLib)

set_property(TARGET
MathVectorTest
MathMatrixTest
Expand All @@ -85,6 +87,8 @@ set_property(TARGET

MathDistanceTest
MathIntersectionTest

MathRandomTest
APPEND PROPERTY COMPILE_DEFINITIONS "CORRADE_GRACEFUL_ASSERT")

set_target_properties(
Expand Down Expand Up @@ -130,4 +134,6 @@ set_target_properties(

MathConfigurationValueTest
MathStrictWeakOrderingTest

MathRandomTest
PROPERTIES FOLDER "Magnum/Math/Test")
101 changes: 101 additions & 0 deletions src/Magnum/Math/Test/RandomTest.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,101 @@
#include <Corrade/TestSuite/Tester.h>
#include <Corrade/TestSuite/Compare/Numeric.h>
#include <Corrade/Utility/DebugStl.h>

#include "Magnum/Math/Random.h"

namespace Magnum
{
namespace Math
{

namespace Test
{
namespace
{

struct RandomTest : Corrade::TestSuite::Tester
{
explicit RandomTest();

void randScalar();
void unitVector2();
void unitVector3();
void randomRotation();
void randomDiceChiSquare();
};

typedef Vector<2, Float> Vector2;
typedef Vector<3, Float> Vector3;
typedef Math::Constants<Float> Constants;

RandomTest::RandomTest()
{
Corrade::TestSuite::Tester::addRepeatedTests(
{&RandomTest::randScalar,
&RandomTest::unitVector2,
&RandomTest::unitVector3,
&RandomTest::randomRotation},
/*repeat number*/ 200);
Corrade::TestSuite::Tester::addTests(
{&RandomTest::randomDiceChiSquare});
}

void RandomTest::randScalar()
{
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()
{
Math::Random::RandomGenerator g;
CORRADE_COMPARE((Math::Random::randomUnitVector2(g)).length(), 1.0f);
}
void RandomTest::unitVector3()
{
Math::Random::RandomGenerator g;

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

void RandomTest::randomRotation()
{
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;
const Int expected = 10000;
const Float thresholdfor100 = 36.191;

for (auto i = 0; i < 100; i++)
{
std::vector<Int> faces(dice_side, 0);
for (std::size_t i = 0; i < expected * dice_side; i++)
faces[Math::Random::randomScalar<Int>(g, 0, dice_side - 1)]++;
Float chi_square = 0.0f;
for (std::size_t i = 0; i < dice_side; i++)
chi_square += Float(pow((faces[i] - expected), 2)) / expected;
if (chi_square > thresholdfor100)
error_count++;
}
CORRADE_COMPARE_AS(error_count, 2, Corrade::TestSuite::Compare::Less);
}
Copy link
Owner

Choose a reason for hiding this comment

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

It would be great to have some sort of distribution verification in the tests, to ensure we're not accidentally skewing the distribution to something non-uniform -- how good are your statistics skills? :)

Found this on SO, it suggests using a Chi Squared test. I never did such a thing myself tho :D

Copy link
Author

Choose a reason for hiding this comment

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

We have easter incoming ! :)

Copy link
Author

Choose a reason for hiding this comment

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

Well, I put a chi square test. In testing, I couldnt see a "tolerant" testing scheme(Chi square let 1 fail; like 99% ).


} // namespace
} // namespace Test
} // namespace Math
} // namespace Magnum

CORRADE_TEST_MAIN(Magnum::Math::Test::RandomTest)