Skip to content

prevent random distributions from returning NaN/Inf #1228

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.

Sign up for GitHub

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

Merged
merged 12 commits into from
Oct 3, 2020
Merged

prevent random distributions from returning NaN/Inf #1228

Changes from 3 commits
Commits
Show all changes
12 commits
Select commit Hold shift + click to select a range
6f7eb15
Change the base to the most recent commit on master
Arzaghi Aug 31, 2020
a0657ed
find the exact lower bound causes overflow
Arzaghi Sep 2, 2020
01ad653
use scaling method to prevent overflow
Arzaghi Sep 7, 2020
207f232
Update stl/inc/random
Arzaghi Sep 9, 2020
89aed7e
Apply suggestions from code review
Arzaghi Sep 9, 2020
da5ca51
Apply suggestions from code review
Arzaghi Sep 9, 2020
91b5ca0
prevent student_t_distribution from returning Inf or NaN
Arzaghi Sep 9, 2020
109ec2a
add cover tests to validate random distributions output for Inf and NaN
Arzaghi Sep 9, 2020
2b124c3
Apply suggestions from code review
Arzaghi Sep 9, 2020
93a68e0
insert an explicit cast
Arzaghi Sep 9, 2020
58c258c
Code review feedback.
StephanTLavavej Oct 2, 2020
7298946
Merge branch 'master' into gh1228
StephanTLavavej Oct 2, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
53 changes: 43 additions & 10 deletions stl/inc/random
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2932,15 +2932,32 @@ private:
_Ty _Vx1;
_Ty _Vx2;
_Ty _Sx;
for (;;) { // reject bad values
for (;;) { // reject bad values to avoid generating NaN/Inf on the next calculations
_Vx1 = 2 * _NRAND(_Eng, _Ty) - 1;
_Vx2 = 2 * _NRAND(_Eng, _Ty) - 1;
_Sx = _Vx1 * _Vx1 + _Vx2 * _Vx2;
if (_Sx < 1) {
if (_Sx < _Ty{1} && _Sx != _Ty{0}) {
// good values!
break;
}
}
const auto _Fx = static_cast<_Ty>(_CSTD sqrt(-2.0 * _CSTD log(_Sx) / _Sx));

_Ty _LogSx;
if (_Sx > 1.0e-4) {
_LogSx = _STD log(_Sx);
} else {
// Bad _Sx value! It will cause overflow on log calculation
// generating a new value based on scaling method.
const _Ty _ln2{static_cast<_Ty>(0.69314718055994530941723212145818)};
const _Ty _Maxabs{_STD max(_STD abs(_Vx1), _STD abs(_Vx2))};
const int expMax{_STD ilogb(_Maxabs)};
_Vx1 = _STD scalbn(_Vx1, -expMax);
_Vx2 = _STD scalbn(_Vx2, -expMax);
_Sx = _Vx1 * _Vx1 + _Vx2 * _Vx2;
_LogSx = _STD log(_Sx) + static_cast<_Ty>(expMax) * (_ln2 * 2);
}

const auto _Fx = static_cast<_Ty>(_STD sqrt(static_cast<_Ty>(-2) * _LogSx / _Sx));
if (_Keep) { // save second value for next call
_Xx2 = _Fx * _Vx2;
_Valid = true;
Expand Down Expand Up @@ -3864,19 +3881,29 @@ public:
_Px1 = _CSTD pow(_Px1, _Ty{1} / _Ax);
_Px2 = _CSTD pow(_Px2, _Ty{1} / _Bx);
_Wx = _Px1 + _Px2;
if (_Wx <= _Ty{1}) {
if (_Wx <= _Ty{1} && _Wx != _Ty{0}) {
break;
}
}
return _Px1 / _Wx;
} else { // use gamma distributions instead
_Ty _Px1;
_Ty _Px2;
_Ty _PSum;
gamma_distribution<_Ty> _Dist1(_Ax, 1);
gamma_distribution<_Ty> _Dist2(_Bx, 1);
_Px1 = _Dist1(_Eng);
_Px2 = _Dist2(_Eng);
return _Px1 / (_Px1 + _Px2);

for (;;) { // reject pairs whose sum is zero
_Px1 = _Dist1(_Eng);
_Px2 = _Dist2(_Eng);
_PSum = _Px1 + _Px2;

if (_PSum != _Ty{0}) {
break;
}
}

return _Px1 / _PSum;
}
}

Expand Down Expand Up @@ -4001,12 +4028,18 @@ private:
_Ty _Px;
_Ty _Vx1;
_Ty _Vx2;
const _Ty _Vx3{1};
_Vx1 = static_cast<_Ty>(_Par0._Mx) * static_cast<_Ty>(0.5);
_Vx2 = static_cast<_Ty>(_Par0._Nx) * static_cast<_Ty>(0.5);
_Beta_distribution<_Ty> _Dist(_Vx1, _Vx2);
_Px = _Dist(_Eng);
for (;;) { // reject bad values
_Px = _Dist(_Eng);
if (_Px != _Vx3) {
break;
}
}

return (_Vx2 / _Vx1) * (_Px / (_Ty{1} - _Px));
return (_Vx2 / _Vx1) * (_Px / (_Vx3 - _Px));
}

param_type _Par;
Expand Down Expand Up @@ -4137,7 +4170,7 @@ private:
_Vx2 = _Dist(_Eng);
_Rs = _Vx1 * _Vx1 + _Vx2 * _Vx2;

if (_Rs < _Ty{1}) {
if (_Rs < _Ty{1} && _Rs != _Ty{0}) {
break;
}
}
Expand Down