diff --git a/python/damask/util.py b/python/damask/util.py index b4ee64283..b5fc2a9d7 100644 --- a/python/damask/util.py +++ b/python/damask/util.py @@ -265,7 +265,7 @@ def show_progress(iterable: _Iterable, N_iter : int, optional Total number of iterations. Required if iterable is not a sequence. prefix : str, optional - Prefix string. + Prefix string. Defaults to '' bar_length : int, optional Length of progress bar in characters. Defaults to 50. @@ -291,7 +291,8 @@ def show_progress(iterable: _Iterable, status.update(i) -def scale_to_coprime(v: _FloatSequence) -> _np.ndarray: +def scale_to_coprime(v: _FloatSequence, + N_significant: int = 9) -> _np.ndarray: """ Scale vector to co-prime (relatively prime) integers. @@ -299,6 +300,8 @@ def scale_to_coprime(v: _FloatSequence) -> _np.ndarray: ---------- v : sequence of float, len (:) Vector to scale. + N_significant: int, optional + Number of significant digits to consider. Defaults to 9 Returns ------- @@ -306,26 +309,25 @@ def scale_to_coprime(v: _FloatSequence) -> _np.ndarray: Vector scaled to co-prime numbers. """ - MAX_DENOMINATOR = 1000000 - def get_square_denominator(x): + def get_square_denominator(x,max_denominator): """Denominator of the square of a number.""" - return _fractions.Fraction(x ** 2).limit_denominator(MAX_DENOMINATOR).denominator + return _fractions.Fraction(x ** 2).limit_denominator(max_denominator).denominator def lcm(a,b): """Least common multiple.""" try: - return _np.lcm(a,b) # numpy > 1.18 + return _np.abs(_np.lcm(a,b)) # numpy > 1.18 except AttributeError: - return a * b // _np.gcd(a, b) + return _np.abs(a * b // _np.gcd(a, b)) - v_ = _np.array(v) - m = (v_ * _reduce(lcm, map(lambda x: int(get_square_denominator(x)),v_))**0.5).astype(_np.int64) + v_ = _np.round(_np.array(v,'float64')/_np.max(_np.abs(v)),N_significant) + max_denominator = int(10**(N_significant-1)) + m = (v_ * _reduce(lcm, map(lambda x: int(get_square_denominator(x,max_denominator)),v_))**0.5).astype(_np.int64) m = m//_reduce(_np.gcd,m) - with _np.errstate(invalid='ignore'): - if not _np.allclose(_np.ma.masked_invalid(v_/m),v_[_np.argmax(abs(v_))]/m[_np.argmax(abs(v_))]): - raise ValueError(f'invalid result "{m}" for input "{v_}"') + if not _np.allclose(m/_np.max(_np.abs(m)),v/_np.max(_np.abs(v)),atol=1e-2,rtol=0): + raise ValueError(f'invalid result "{m}" for input "{v}"') return m diff --git a/python/tests/test_util.py b/python/tests/test_util.py index ce5d1b3a9..0eac5febb 100644 --- a/python/tests/test_util.py +++ b/python/tests/test_util.py @@ -37,20 +37,20 @@ class TestUtil: def test_srepr(self,input,glue,quote,output): assert output == util.srepr(input,glue,quote) + + @pytest.mark.parametrize('N',[5,6,7,8,9,10,11,12,13,14,15,16,20,30,40,50]) @pytest.mark.parametrize('input,output', [ ([0,-2],[0,-1]), ([-0.5,0.5],[-1,1]), ([1./2.,1./3.],[3,2]), ([2./3.,1./2.,1./3.],[4,3,2]), + ([0.666666666666,-0.33333333333,-0.33333],[2,-1,-1]), + ([1./3., 1./4., 1./22],[536870912, 402653184, 73209669]), ]) - def test_scale2coprime(self,input,output): - assert np.allclose(util.scale_to_coprime(np.array(input)), - np.array(output).astype(int)) - - def test_lackofprecision(self): - with pytest.raises(ValueError): - util.scale_to_coprime(np.array([1/333.333,1,1])) + def test_scale2coprime(self,input,output,N): + res = util.scale_to_coprime(input,N) + assert np.allclose(res/np.max(np.abs(res)),output/np.max(np.abs(output)),atol=1e-2,rtol=0) @pytest.mark.parametrize('rv',[stats.rayleigh(),stats.weibull_min(1.2),stats.halfnorm(),stats.pareto(2.62)])