Coverage for pygeodesy/fmath.py : 85%

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 -*-
@newfield example: Example, Examples ''' # make sure int/int division yields float quotient raise ImportError('%s 1/2 == %d' % ('division', division))
# all public contants, classes and functions
except ImportError: try: # _Ints imported by .utily _Ints = int, long #: (INTERNAL) Int objects (C{tuple}) except NameError: # Python 3+ _Ints = int, #: (INTERNAL) Int objects (C{tuple})
except ImportError: try: _Scalars = int, long, float #: (INTERNAL) Scalar objects (C{tuple}) except NameError: _Scalars = int, float #: (INTERNAL) Scalar objects (C{tuple})
except ImportError: _Seqs = list, tuple, range # XXX also set?
except AttributeError: EPS = 2.220446049250313e-16 #: Epsilon (C{float}) 2**-52? MANTIS = 53 #: Mantissa bits ≈53 (C{int}) MAX = pow(2.0, 1023) * (2 - EPS) #: Float max (C{float}) ≈10**308, 2**1024? MIN = pow(2.0, -1021) # Float min (C{float}) ≈10**-308, 2**-1021? # _1EPS = 1.0 + EPS #: M{1 + EPS} ≈1.0000000000000002 (C{float})
'''(INTERNAL) Format a C{TypeError} for a C{name=value} pair. ''' else: n, v = 'pair', 'N/A'
'''(INTERNAL) Half-even rounding. ''' (r < 0 and p < 0): # signs match
'''(INTERNAL) Precision C{2sum} of M{a + b}. ''' raise OverflowError('%s: %r' % (_2sum.__name__, s))
'''Precision summation similar to standard Python function C{math.fsum}.
Unlike C{math.fsum}, this class accumulates the values repeatedly and provides intermediate, precision running sums. Accumulation may continue after intermediate summations.
@note: Handling of exceptions, C{nan} and C{finite} values is different from C{math.fsum}.
@see: U{Hettinger<https://code.ActiveState.com/recipes/393090>}, 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.
@param 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.
@param 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.
@param 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('%s += %r' % (self, other))
'''Multiply this instance by a scalar or an other instance.
@param 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}. ''' p = s.fcopy() p.fmul(ps.pop()) self.fadd(p._ps) else: self._ps = [] # zero self._fsum2_ = None else: raise TypeError('%s *= %r' % (self, other))
'''Subtract a scalar or an other instance from this instance.
@param 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('%s -= %r' % (self, other))
'''Return the number of accumulated values. '''
'''Product of this and an other instance or a scalar.
@param other: L{Fsum} instance or C{scalar}.
@return: The product, a new instance (L{Fsum}).
@see: Method L{Fsum.__imul__}. '''
# m = self.__module__.split('.')[-1]
'''Difference of this and an other instance or a scalar.
@param 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.
@param iterable: Sequence, list, tuple, etc. (C{scalar}s).
@raise OverflowError: Partial C{2sum} overflow.
@raise TypeError: Non-scalar B{C{iterable}} value.
@raise ValueError: Invalid or non-finite B{C{iterable}} value. ''' iterable = tuple(iterable)
# def _iter(): # for a in iterable: # if isinstance(a, Fsum): # if a is self: # self.fmul(2) # else: # for a in a._ps: # yield a # else: # yield a
raise ValueError('%s, not %s: %r' % (self, 'finite', a)) # assert self._ps is ps
'''Accumulate more values from positional arguments.
@param xs: Values to add (C{scalar}s), all positional.
@see: Method L{Fsum.fadd}. '''
'''Copy this instance, shallow or deep.
@return: The copy, a new instance (L{Fsum}). '''
'''Multiple the current, partial sum by a factor.
@param 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('%s, not %s: %r' % (self, 'finite', factor))
# assert self._ps is ps
'''Accumulate more values from an iterable.
@param iterable: Sequence, list, tuple, etc. (C{scalar}s).
@see: Method L{Fsum.fadd}. '''
'''Accumulate more values from positional arguments.
@param xs: Values to subtract (C{scalar}s), all positional.
@see: Method L{Fsum.fadd}. '''
'''Accumulate more values from an iterable and sum all.
@keyword 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.
@param 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.
@param 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))}.
@param a: List, sequence, tuple, etc. (C{scalar}s). @param b: All positional arguments (C{scalar}s).
@raise OverflowError: Partial C{2sum} overflow.
@raise ValueError: Unequal C{len}(B{a}) and C{len}(B{b}).
@see: Function L{fdot} and method L{Fsum.fadd}. ''' raise ValueError('%s, %s: %s vs %s' % (self, 'len', len(a), 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))}.
@param x: Polynomial argument (C{scalar}). @param 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('%s, not %s: %r' % (self, 'finite', x)) raise ValueError('%s, no %s: %r' % (self, 'coefficents', cs))
'''Precision polynomial evaluation. ''' '''New L{Fpolynomial} evaluation of the polynomial M{sum(cs[i] * x**i for i=0..len(cs))}.
@param x: Polynomial argument (C{scalar}). @param 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('%s, not %s: %r' % (self, 'finite', x)) raise ValueError('%s, no %s: %r' % (self, 'coefficents', cs))
'''Return M{math.acos(max(-1, min(1, x)))}. '''
'''Compute the cubic root M{x**(1/3)}.
@param 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)}.
@param x: Value (C{scalar}).
@return: Cubic root squared (C{float}).
@see: Functions L{cbrt} and L{sqrt3}. '''
'''Return the weighted average of two values.
@param v1: One value (C{scalar}). @param v2: Other value (C{scalar}). @keyword 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('%s invalid: %r' % ('fraction', f))
'''Return the precision dot product M{sum(a[i] * b[i] for i=0..len(a))}.
@param a: List, sequence, tuple, etc. (C{scalar}s). @param b: All positional arguments (C{scalar}s).
@return: Dot product (C{float}).
@raise ValueError: Unequal C{len(B{a})} and C{len(B{b})}.
@see: Class L{Fdot}. ''' raise ValueError('%s(%s): %s vs %s' % (fdot.__name__, 'len', len(a), len(b)))
'''Return the precision dot product M{start + sum(a[i] * b[i] * c[i] for i=0..len(a))}.
@param a: List, sequence, tuple, etc. (C{scalar}[]). @param b: List, sequence, tuple, etc. (C{scalar}[]). @param c: List, sequence, tuple, etc. (C{scalar}[]). @keyword start: Optional bias (C{scalar}).
@return: Dot product (C{float}).
@raise OverflowError: Partial C{2sum} overflow.
@raise ValueError: Unequal C{len(B{a})}, C{len(B{b})} and/or C{len(B{c})}. '''
raise ValueError('%s(%s): %s vs %s vs %s' % (fdot3.__name__, 'len', len(a), len(b), len(c)))
else:
'''Evaluate the polynomial M{sum(cs[i] * x**i for i=0..len(cs))} using the Horner form.
@param x: Polynomial argument (C{scalar}). @param 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).
@param xs: Known values (C{scalar}[]). @param ds: Non-negative distances (C{scalar}[]). @keyword beta: Inverse distance power (C{int}, 0, 1, 2, or 3).
@return: Interpolated value C{x} (C{float}).
@raise ValueError: Invalid B{C{beta}}, negative B{C{ds}} value, weighted B{C{ds}} below L{EPS} or unequal C{len(B{ds})} and C{len(B{xs})}.
@note: Using C{B{beta}=0} returns the mean of B{C{xs}}. ''' raise ValueError('%s(%s): %s vs %s' % (fidw.__name, 'len', n, d))
raise ValueError('%s(%s[%s]) invalid: %r' % (fidw.__name, 'ds', '', d)) elif b == 0: x = fmean(xs) else: raise ValueError('%s(%s=%r) invalid' % (fidw.__name, 'beta', beta)) i = ds.index(d) raise ValueError('%s(%s[%s]) invalid: %r' % (fidw.__name, 'ds', i, d))
'''Compute the accurate mean M{sum(xs[i] for i=0..len(xs)) / len(xs)}.
@param 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('%s(%r)' % (fmean.__name__, xs))
'''Evaluate the polynomial M{sum(cs[i] * x**i for i=0..len(cs))}.
@param x: Polynomial argument (C{scalar}). @param 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]}.
@param x: Value (C{scalar}). @param n: Highest exponent (C{int}). @keyword alts: Only alternating powers, starting with this exponent (C{int}).
@return: Powers of B{C{x}} (C{float}[]).
@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('not %s: %r' %('finite', x)) raise _IsNotError(int.__name_, n=n) raise ValueError('%s invalid: %r' % ('n', n))
# XXX PyChecker chokes on xs[alts-1::2]
# XXX PyChecker claims result is None
'''Iterable product, like C{math.prod} or C{numpy.prod}.
@param iterable: Values to be multiplied (C{scalar}[]). @keyword start: Initial product, 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.
@param start: First value (C{float}). @param number: The number of C{float}s to generate (C{int}). @keyword 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: try: freduce = reduce # PYCHOK expected except NameError: # Python 3+ _EMPTY = object()
def freduce(f, iterable, *start): '''For missing C{functools.reduce}. ''' if start: r = v = start[0] else: r, v = 0, _EMPTY for v in iterable: r = f(r, v) if v is _EMPTY: raise TypeError('%s() empty, no start' % (freduce.__name__,)) return r
'''Convert floats to string, optionally with trailing zero decimals stripped.
@param floats: List, sequence, tuple, etc. (C{scalar}s). @keyword prec: Optional precision, number of decimal digits (0..9). Trailing zero decimals are stripped for B{C{prec}} values of 1 and above, but kept for negative B{C{prec}} values. @keyword sep: Optional, separator to join (string). @keyword fmt: Optional, float format (string). @keyword ints: Optionally, remove decimal dot (C{bool}).
@return: The floats as 'f, f, ... f' (string). ''' # corner case testLcc lon0=-96.0 t.rstrip('0').endswith('.')):
else:
'''Strip trailing zero decimals from a float string.
@param fstr: Float (string).
@return: Float (string). '''
'''Precision summation of the positional argument vulues.
@param 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 # no, remove fsum ... raise ImportError # ... use fsum below
except ImportError:
def fsum(iterable): '''Precision summation similar to standard Python function C{math.fsum}.
Exception and I{non-finite} handling differs from C{math.fsum}.
@param 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)
'''Compute the norm M{sqrt(1 + x**2)}.
@param x: Argument (C{scalar}).
@return: Norm (C{float}). '''
'''Compute the norm M{sqrt(sum(xs[i]**2)) for i=0..len(xs)}.
@param 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+ 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+ U{math.dist <https://docs.Python.org/3.8/library/math.html#math.dist>} Euclidian distance 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}. ''' else: n = '' raise ValueError('%s(): %r[%s]' % (hypot_.__name__, xs, n))
except ImportError:
def isfinite(obj): '''Check for C{Inf} and C{NaN} values.
@param obj: Value (C{scalar}).
@return: C{False} if B{C{obj}} is C{INF} or C{NAN}, C{True} otherwise.
@raise TypeError: Non-scalar B{C{obj}}. ''' if not isscalar(obj): raise _IsNotError(isscalar.__name__, obj=obj) return not (isinf(obj) or isnan(obj))
'''Check for integer type or integer value.
@param obj: The object (any C{type}). @keyword both: Optionally, check both type and value (C{bool}).
@return: C{True} if B{C{obj}} is C{int}, C{False} otherwise. ''' except AttributeError: return False # XXX float(int(obj)) == obj?
'''Check for NEG0, negative 0.0.
@param obj: Value (C{scalar}).
@return: C{True} if B{C{obj}} is C{NEG0} or -0.0, C{False} otherwise. ''' # and str(obj).rstrip('0') == '-0.'
'''Check for scalar types.
@param obj: The object (any C{type}).
@return: C{True} if B{C{obj}} is C{scalar}, C{False} otherwise. '''
'''Make built-in function L{len} work for generators, iterators, etc. since those can only be started exactly once.
@param seq: Generator, iterator, list, range, tuple, etc.
@return: 2-Tuple (number, ...) of items (C{int}, C{list} or C{range} or C{tuple}). '''
'''Apply each argument to a single-argument function and return a tuple of results.
@param func: Function to apply (C{callable}). @param xs: Arguments to apply (C{any positional}).
@return: Function results (C{tuple}). '''
'''Apply arguments to a function and return a tuple of results.
Unlike Python 2's built-in L{map}, Python 3+ L{map} returns a L{map} object, an iterator-like object which generates the results only once. Converting the L{map} object to a tuple maintains Python 2 behavior.
@param func: Function to apply (C{callable}). @param xs: Arguments to apply (C{list, tuple, ...}).
@return: Function results (C{tuple}). '''
'''Validate a scalar.
@param value: The value (C{scalar}). @keyword low: Optional lower bound (C{scalar}). @keyword high: Optional upper bound (C{scalar}). @keyword name: Optional name of value (C{str}). @keyword Error: Exception to raise (C{ValueError}).
@return: New value (C{type} of B{C{low}}).
@raise TypeError: Non-scalar B{C{value}}.
@raise Error: Out-of-bounds B{C{value}}. ''' raise _IsNotError(scalar.__name__, **{name: value}) else: raise ValueError except (TypeError, ValueError): raise _IsNotError('valid', Error=Error, **{name: value})
'''Compute the square root cubed M{sqrt(x)**3} or M{sqrt(x**3)}.
@param 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('%s(%r)' % (sqrt3.__name__, x))
# **) MIT License # # Copyright (C) 2016-2020 -- 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. |