Coverage for pygeodesy/fmath.py : 95%

Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
# -*- coding: utf-8 -*-
''' # make sure int/int division yields float quotient, see .basics
isscalar, len2, _xcopy _TypeError, _ValueError _EPS0__2, _finite_, _few_, _h_, _negative_,\ _not_, _singular_, _SPACE_, _too_, \ _0_0, _1_0, _1_5 as _3_2nd, _2_0, _3_0
# sqrt(2) <https://WikiPedia.org/wiki/Square_root_of_2>
'''(INTERNAL) Half-even rounding. ''' (r < 0 and p < 0): # signs match
'''(INTERNAL) Precision C{2sum} of M{a + b}. ''' raise _OverflowError(unstr(_2sum.__name__, a, b), txt=str(s))
'''Precision summation similar to standard Python function C{math.fsum}.
Unlike C{math.fsum}, this class accumulates the values and provides intermediate, precision running sums. Accumulation may continue after intermediate summations.
@note: Handling of exceptions, C{inf}, C{INF}, C{nan} and C{NAN} values is different from C{math.fsum}.
@see: U{Hettinger<https://GitHub.com/ActiveState/code/blob/master/recipes/Python/ 393090_Binary_floating_point_summatiaccurate_full/recipe-393090.py>}, U{Kahan<https://WikiPedia.org/wiki/Kahan_summation_algorithm>}, U{Klein<https://Link.Springer.com/article/10.1007/s00607-005-0139-x>}, Python 2.6+ file I{Modules/mathmodule.c} and the issue log U{Full precision summation<https://Bugs.Python.org/issue2819>}. '''
'''Initialize a new accumulator with one or more start values.
@arg starts: No, one or more start values (C{scalar}s).
@raise OverflowError: Partial C{2sum} overflow.
@raise TypeError: Non-scalar B{C{starts}} value.
@raise ValueError: Invalid or non-finite B{C{starts}} value. '''
'''Sum of this and an other instance or a scalar.
@arg other: L{Fsum} instance or C{scalar}.
@return: The sum, a new instance (L{Fsum}).
@see: Method L{Fsum.__iadd__}. '''
'''Add a scalar or an other instance to this instance.
@arg other: L{Fsum} instance or C{scalar}.
@return: This instance, updated (L{Fsum}).
@raise TypeError: Invalid B{C{other}} type.
@see: Method L{Fsum.fadd}. ''' else: raise _TypeError(_SPACE_(self, '+=', repr(other)))
'''Multiply this instance by a scalar or an other instance.
@arg other: L{Fsum} instance or C{scalar}.
@return: This instance, updated (L{Fsum}).
@raise TypeError: Invalid B{C{other}} type.
@see: Method L{Fsum.fmul}. ''' else: self._ps = [] # zero self._fsum2_ = None else: raise _TypeError(_SPACE_(self, '*=', repr(other)))
'''Subtract a scalar or an other instance from this instance.
@arg other: L{Fsum} instance or C{scalar}.
@return: This instance, updated (L{Fsum}).
@raise TypeError: Invalid B{C{other}} type.
@see: Method L{Fsum.fadd}. ''' else: raise _TypeError(_SPACE_(self, '-=', repr(other)))
'''Return the number of accumulated values. '''
'''Product of this and an other instance or a scalar.
@arg other: L{Fsum} instance or C{scalar}.
@return: The product, a new instance (L{Fsum}).
@see: Method L{Fsum.__imul__}. '''
'''Difference of this and an other instance or a scalar.
@arg other: L{Fsum} instance or C{scalar}.
@return: The difference, a new instance (L{Fsum}).
@see: Method L{Fsum.__isub__}. '''
'''Accumulate more values from an iterable.
@arg xs: Iterable, list, tuple, etc. (C{scalar}s).
@raise OverflowError: Partial C{2sum} overflow.
@raise TypeError: Non-scalar B{C{xs}} value.
@raise ValueError: Invalid or non-finite B{C{xs}} value. ''' xs = (xs,)
n = Fmt.SQUARE(xs=n) raise _ValueError(n, x, txt=_not_(_finite_)) # assert self._ps is ps
'''Accumulate more values from positional arguments.
@arg xs: Values to add (C{scalar}s), all positional.
@raise OverflowError: Partial C{2sum} overflow.
@raise TypeError: Non-scalar B{C{xs}} value.
@raise ValueError: Invalid or non-finite B{C{xs}} value. '''
'''Copy this instance, C{shallow} or B{C{deep}}.
@return: The copy, a new instance (L{Fsum}). ''' # f._fsum2_ = self._fsum2_ # f._n = self._n
'''Multiple the current, partial sum by a factor.
@arg factor: The multiplier (C{scalar}).
@raise TypeError: Non-scalar B{C{factor}}.
@raise ValueError: Invalid or non-finite B{C{factor}}.
@see: Method L{Fsum.fadd}. ''' raise _ValueError(factor=factor, txt=_not_(_finite_))
# assert self._ps is ps
'''Accumulate more values from an iterable.
@arg iterable: Sequence, list, tuple, etc. (C{scalar}s).
@see: Method L{Fsum.fadd}. '''
'''Accumulate more values from positional arguments.
@arg xs: Values to subtract (C{scalar}s), all positional.
@see: Method L{Fsum.fadd}. '''
'''Accumulate more values from an iterable and sum all.
@kwarg iterable: Sequence, list, tuple, etc. (C{scalar}s), optional.
@return: Accurate, running sum (C{float}).
@raise OverflowError: Partial C{2sum} overflow.
@raise TypeError: Non-scalar B{C{iterable}} value.
@raise ValueError: Invalid or non-finite B{C{iterable}} value.
@note: Accumulation can continue after summation. '''
else: # assert self._ps is ps
'''Accumulate more values from positional arguments and sum all.
@arg xs: Values to add (C{scalar}s), all positional.
@return: Accurate, running sum (C{float}).
@see: Method L{Fsum.fsum}.
@note: Accumulation can continue after summation. '''
'''Accumulate more values from positional arguments, sum all and provide the sum and delta.
@arg xs: Values to add (C{scalar}s), all positional.
@return: 2-Tuple C{(sum, delta)} with the accurate, running C{sum} and the C{delta} with the previous running C{sum}, both (C{float}).
@see: Method L{Fsum.fsum_}.
@note: Accumulation can continue after summation. '''
'''Precision dot product. ''' '''New L{Fdot} precision dot product M{sum(a[i] * b[i] for i=0..len(a))}.
@arg a: List, sequence, tuple, etc. (C{scalar}s). @arg b: All positional arguments (C{scalar}s).
@raise OverflowError: Partial C{2sum} overflow.
@raise LenError: Unequal C{len(B{a})} and C{len(B{b})}.
@see: Function L{fdot} and method L{Fsum.fadd}. ''' raise LenError(Fdot, a=len(a), b=len(b))
'''Precision polynomial evaluation using the Horner form. ''' '''New L{Fhorner} evaluation of the polynomial M{sum(cs[i] * x**i for i=0..len(cs))}.
@arg x: Polynomial argument (C{scalar}). @arg cs: Polynomial coeffients (C{scalar}[]).
@raise OverflowError: Partial C{2sum} overflow.
@raise TypeError: Non-scalar B{C{x}}.
@raise ValueError: No B{C{cs}} coefficients or B{C{x}} is not finite.
@see: Function L{fhorner} and methods L{Fsum.fadd} and L{Fsum.fmul}. ''' raise _ValueError(x=x, txt=_not_(_finite_)) raise _ValueError(cs=cs, txt=MISSING)
'''Precision polynomial evaluation. ''' '''New L{Fpolynomial} evaluation of the polynomial M{sum(cs[i] * x**i for i=0..len(cs))}.
@arg x: Polynomial argument (C{scalar}). @arg cs: Polynomial coeffients (C{scalar}[]).
@raise OverflowError: Partial C{2sum} overflow.
@raise TypeError: Non-scalar B{C{x}}.
@raise ValueError: No B{C{cs}} coefficients or B{C{x}} is not finite.
@see: Function L{fpolynomial} and method L{Fsum.fadd}. ''' raise _ValueError(x=x, txt=_not_(_finite_)) raise _ValueError(cs=cs, txt=MISSING)
'''Compute the cubic root M{x**(1/3)}.
@arg x: Value (C{scalar}).
@return: Cubic root (C{float}).
@see: Functions L{cbrt2} and L{sqrt3}. ''' # simpler and more accurate than Ken Turkowski's CubeRoot, see # <https://People.FreeBSD.org/~lstewart/references/apple_tr_kt32_cuberoot.pdf>
'''Compute the cubic root squared M{x**(2/3)}.
@arg x: Value (C{scalar}).
@return: Cubic root squared (C{float}).
@see: Functions L{cbrt} and L{sqrt3}. '''
'''I{Appoximate} the norm M{sqrt(x**2 + y**2)} by M{max(abs(x), abs(y)) + min(abs(x), abs(y)) * 0.4142...}.
@arg x: X component (C{scalar}). @arg y: Y component (C{scalar}).
@return: Appoximate norm (C{float}).
@see: Function L{euclid_}. '''
'''I{Appoximate} the norm M{sqrt(sum(x**2 for x in xs))} by cascaded L{euclid}.
@arg xs: X arguments, positional (C{scalar}[]).
@return: Appoximate norm (C{float}).
@see: Function L{euclid}. ''' # e = euclid(x, e)
'''Fast approximation of L{acos1}C{(B{x})}.
@see: U{ShaderFastLibs.h<https://GitHub.com/michaldrobot/ ShaderFastLibs/blob/master/ShaderFastMathLib.h>}. ''' r = PI_2 else:
def fasin1(x): # PYCHOK no cover '''Fast approximation of L{asin1}C{(B{x})}.
@see: L{facos1}. ''' return PI_2 - facos1(x)
'''Fast approximation of C{atan(B{x})}. ''' elif a > _1_0: r = PI_2 - fatan1(_1_0 / a) # == fatan2(a, _1_0) else: r = PI_4
'''Fast approximation of C{atan(B{x})} for C{0 <= B{x} <= 1}, I{unchecked}.
@see: U{ShaderFastLibs.h<https://GitHub.com/michaldrobot/ ShaderFastLibs/blob/master/ShaderFastMathLib.h>} and U{Efficient approximations for the arctangent function <http://www-Labs.IRO.UMontreal.CA/~mignotte/IFT2425/Documents/ EfficientApproximationArctgFunction.pdf>}, IEEE Signal Processing Magazine, 111, May 2006. ''' # Eq (9): PI_4 * x - x * (x - 1) * (0.2447 + 0.0663 * x**2)
'''Fast approximation of C{atan2(B{y}, B{x})}.
@see: U{fastApproximateAtan(x, y)<https://GitHub.com/CesiumGS/cesium/blob/ master/Source/Shaders/Builtin/Functions/fastApproximateAtan.glsl>} and L{fatan1}. ''' else: # a == b == 0 return _0_0 r = PI - r
'''Return the average of two values.
@arg v1: One value (C{scalar}). @arg v2: Other value (C{scalar}). @kwarg f: Optional fraction (C{float}).
@return: M{v1 + f * (v2 - v1)} (C{float}). ''' # @raise ValueError: Fraction out of range. # ''' # if not 0 <= f <= 1: # XXX restrict fraction? # raise _ValueError(fraction=f)
'''Return the precision dot product M{sum(a[i] * b[i] for i=0..len(a))}.
@arg a: List, sequence, tuple, etc. (C{scalar}s). @arg b: All positional arguments (C{scalar}s).
@return: Dot product (C{float}).
@raise LenError: Unequal C{len(B{a})} and C{len(B{b})}.
@see: Class L{Fdot}. ''' raise LenError(fdot, a=len(a), b=len(b))
'''Return the precision dot product M{start + sum(a[i] * b[i] * c[i] for i=0..len(a))}.
@arg a: List, sequence, tuple, etc. (C{scalar}[]). @arg b: List, sequence, tuple, etc. (C{scalar}[]). @arg c: List, sequence, tuple, etc. (C{scalar}[]). @kwarg start: Optional bias (C{scalar}).
@return: Dot product (C{float}).
@raise LenError: Unequal C{len(B{a})}, C{len(B{b})} and/or C{len(B{c})}.
@raise OverflowError: Partial C{2sum} overflow. '''
raise LenError(fdot3, a=len(a), b=len(b), c=len(c))
else:
'''Evaluate the polynomial M{sum(cs[i] * x**i for i=0..len(cs))} using the Horner form.
@arg x: Polynomial argument (C{scalar}). @arg cs: Polynomial coeffients (C{scalar}[]).
@return: Horner value (C{float}).
@raise OverflowError: Partial C{2sum} overflow.
@raise TypeError: Non-scalar B{C{x}}.
@raise ValueError: No B{C{cs}} coefficients or B{C{x}} is not finite.
@see: Function L{fpolynomial} and class L{Fhorner}. '''
'''Interpolate using using U{Inverse Distance Weighting <https://WikiPedia.org/wiki/Inverse_distance_weighting>} (IDW).
@arg xs: Known values (C{scalar}[]). @arg ds: Non-negative distances (C{scalar}[]). @kwarg beta: Inverse distance power (C{int}, 0, 1, 2, or 3).
@return: Interpolated value C{x} (C{float}).
@raise LenError: Unequal or zero C{len(B{ds})} and C{len(B{xs})}.
@raise ValueError: Invalid B{C{beta}}, negative B{C{ds}} value, weighted B{C{ds}} below L{EPS}.
@note: Using C{B{beta}=0} returns the mean of B{C{xs}}. ''' raise LenError(fidw, xs=n, ds=d)
n = Fmt.PAREN(fsum='ds') raise _ValueError(n, d, txt=_singular_) else: # b == 0 x = fsum(xs) / n # fmean(xs) n = Fmt.SQUARE(ds=ds.index(d)) raise _ValueError(n, d, txt=_negative_)
'''Compute the accurate mean M{sum(xs[i] for i=0..len(xs)) / len(xs)}.
@arg xs: Values (C{scalar}s).
@return: Mean value (C{float}).
@raise OverflowError: Partial C{2sum} overflow.
@raise ValueError: No B{C{xs}} values. ''' raise _ValueError(xs=xs)
'''Compute the accurate mean M{sum(xs[i] for i=0..len(xs)) / len(xs)}.
@arg xs: Values (C{scalar}s).
@return: Mean value (C{float}).
@raise OverflowError: Partial C{2sum} overflow.
@raise ValueError: No B{C{xs}} values. '''
'''Evaluate the polynomial M{sum(cs[i] * x**i for i=0..len(cs))}.
@arg x: Polynomial argument (C{scalar}). @arg cs: Polynomial coeffients (C{scalar}[]).
@return: Polynomial value (C{float}).
@raise OverflowError: Partial C{2sum} overflow.
@raise TypeError: Non-scalar B{C{x}}.
@raise ValueError: No B{C{cs}} coefficients or B{C{x}} is not finite.
@see: Function L{fhorner} and class L{Fpolynomial}. '''
'''Return a series of powers M{[x**i for i=1..n]}.
@arg x: Value (C{scalar}). @arg n: Highest exponent (C{int}). @kwarg alts: Only alternating powers, starting with this exponent (C{int}).
@return: Powers of B{C{x}} (C{float}s or C{int}s).
@raise TypeError: Non-scalar B{C{x}} or B{C{n}} not C{int}.
@raise ValueError: Non-finite B{C{x}} or non-positive B{C{n}}. ''' raise _ValueError(x=x, txt=_not_(_finite_)) raise _IsnotError(int.__name__, n=n) raise _ValueError(n=n)
# xs[alts-1::2] chokes PyChecker
'''Iterable product, like C{math.prod} or C{numpy.prod}.
@arg iterable: Terms to be multiplied (C{scalar}[]). @kwarg start: Initial term, also the value returned for an empty iterable (C{scalar}).
@return: The product (C{float}).
@see: U{NumPy.prod<https://docs.SciPy.org/doc/ numpy/reference/generated/numpy.prod.html>}. ''' return freduce(_mul, iterable, start)
'''Generate a range of C{float}s.
@arg start: First value (C{float}). @arg number: The number of C{float}s to generate (C{int}). @kwarg step: Increment value (C{float}).
@return: A generator (C{float}s).
@see: U{NumPy.prod<https://docs.SciPy.org/doc/ numpy/reference/generated/numpy.arange.html>}. ''' raise _IsnotError(int.__name__, number=number)
except ImportError: # PYCHOK no cover try: freduce = reduce # PYCHOK expected except NameError: # Python 3+
def freduce(f, iterable, *start): '''For missing C{functools.reduce}. ''' if start: r = v = start[0] else: r, v = 0, MISSING for v in iterable: r = f(r, v) if v is MISSING: raise _TypeError(iterable=(), start=MISSING) return r
'''Precision summation of all positional arguments.
@arg xs: Values to be added (C{scalar}[]).
@return: Accurate L{fsum} (C{float}).
@raise OverflowError: Partial C{2sum} overflow.
@raise TypeError: Non-scalar B{C{xs}} value.
@raise ValueError: Invalid or non-finite B{C{xs}} value. '''
# make sure fsum works as expected (XXX check # float.__getformat__('float')[:4] == 'IEEE'?) del fsum # nope, remove fsum ... raise ImportError # ... use fsum below
except ImportError: # PYCHOK no cover
def fsum(iterable): '''Precision summation similar to standard Python function C{math.fsum}.
Exception and I{non-finite} handling differs from C{math.fsum}.
@arg iterable: Values to be added (C{scalar}[]).
@return: Accurate C{sum} (C{float}).
@raise OverflowError: Partial C{2sum} overflow.
@raise TypeError: Non-scalar B{C{iterable}} value.
@raise ValueError: Invalid or non-finite B{C{iterable}} value.
@see: Class L{Fsum}. ''' f = Fsum() return f.fsum(iterable)
from math import hypot # OK in Python 3.10+ hypot_ = hypot else: # Python 3.8 and 3.9 C{math.hypot} is inaccurate, see # agdhruv <https://GitHub.com/geopy/geopy/issues/466>, # cffk <https://Bugs.Python.org/issue43088> and module # geomath.py <https://PyPI.org/project/geographiclib/1.52> else: '''Compute the norm M{sqrt(x**2 + y**2)}.
@arg x: Argument (C{scalar}). @arg y: Argument (C{scalar}).
@return: C{sqrt(B{x}**2 + B{y}**2)} (C{float}). ''' else:
'''Compute the norm M{sqrt(sum(x**2 for x in xs))}.
@arg xs: X arguments, positional (C{scalar}[]).
@return: Norm (C{float}).
@raise OverflowError: Partial C{2sum} overflow.
@raise ValueError: Invalid or no B{C{xs}} value.
@see: Similar to Python 3.8+ n-dimension U{math.hypot <https://docs.Python.org/3.8/library/math.html#math.hypot>}, but handling of exceptions, C{nan} and C{infinite} values is different.
@note: The Python 3.8+ Euclidian distance U{math.dist <https://docs.Python.org/3.8/library/math.html#math.dist>} between 2 I{n}-dimensional points I{p1} and I{p2} can be computed as M{hypot_(*((c1 - c2) for c1, c2 in zip(p1, p2)))}, provided I{p1} and I{p2} have the same, non-zero length I{n}. '''
'''(INTERNAL) Helper for L{hypot_} and L{hypot2_}. '''
else: raise _ValueError(xs=xs, txt=_too_(_few_))
'''Compute the norm M{sqrt(1 + x**2)}.
@arg x: Argument (C{scalar}).
@return: Norm (C{float}). '''
'''Compute the I{squared} norm M{x**2 + y**2}.
@arg x: Argument (C{scalar}). @arg y: Argument (C{scalar}).
@return: C{B{x}**2 + B{y}**2} (C{float}). ''' else:
'''Compute the I{squared} norm C{sum(x**2 for x in B{xs})}.
@arg xs: X arguments, positional (C{scalar}[]).
@return: Squared norm (C{float}).
@raise OverflowError: Partial C{2sum} overflow.
@raise ValueError: Invalid or no B{C{xs}} value.
@see: Function L{hypot_}. '''
'''Normalize a 2-dimensional vector.
@arg x: X component (C{scalar}). @arg y: Y component (C{scalar}).
@return: 2-Tuple C{(x, y)}, normalized.
@raise ValueError: Invalid B{C{x}} or B{C{y}} or zero norm. ''' except (TypeError, ValueError) as X: raise _ValueError(x=x, y=y, h=h, txt=str(X))
'''Normalize all n-dimensional vector components.
@arg xs: The component (C{scalar}[]).
@return: Yield each component, normalized.
@raise ValueError: Invalid or insufficent B{C{xs}} or zero norm. ''' except (TypeError, ValueError) as X: raise _ValueError(Fmt.SQUARE(xs=i), x, _h_, h, txt=str(X))
'''Compute the square root iff C{B{x} > EPS0**2}.
@arg x: Value (C{scalar}).
@return: Square root (C{float}) or C{0.0}. '''
'''Compute the square root, cubed M{sqrt(x)**3} or M{sqrt(x**3)}.
@arg x: Value (C{scalar}).
@return: Cubed square root (C{float}).
@raise ValueError: Negative B{C{x}}.
@see: Functions L{cbrt} and L{cbrt2}. ''' raise _ValueError(x=x)
# **) MIT License # # Copyright (C) 2016-2021 -- mrJean1 at Gmail -- All Rights Reserved. # # Permission is hereby granted, free of charge, to any person obtaining a # copy of this software and associated documentation files (the "Software"), # to deal in the Software without restriction, including without limitation # the rights to use, copy, modify, merge, publish, distribute, sublicense, # and/or sell copies of the Software, and to permit persons to whom the # Software is furnished to do so, subject to the following conditions: # # The above copyright notice and this permission notice shall be included # in all copies or substantial portions of the Software. # # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS # OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL # THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR # OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, # ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR # OTHER DEALINGS IN THE SOFTWARE. |