onyx.util.mathutils – Numerical utilitites

Module attribute ‘primes’ can be used as an (unbounded) generator of prime numbers, or as an (unbounded) sequence of prime numbers.

Look at some primes:

>>> tuple(islice(primes, 10))
(2, 3, 5, 7, 11, 13, 17, 19, 23, 29)
>>> tuple(islice(primes, 20))
(2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71)

The hundredth prime via iteration

>>> tuple(islice(primes, 100))[-1]
541

The thousandth prime via indexing

>>> primes[999]
7919

This module provides some routines for doing log-domain computations on Numpy arrays.

>>> x = np.zeros((3,2), dtype=float) - 3.2
>>> x
array([[-3.2, -3.2],
       [-3.2, -3.2],
       [-3.2, -3.2]])

By default, logsumexp_array does the summation over the entire array, just like Numpy’s sum().

>>> logsumexp_array(x) 
-1.40824053077...

Here’s summation along axis 1, which is the horizontal axis in a 2-d array.

>>> logsumexp_array(x, axis=1)
array([-2.50685282, -2.50685282, -2.50685282])

To do the summation vertically, specify axis 0.

>>> logsumexp_array(x, axis=0)
array([-2.10138771, -2.10138771])

-inf is the additive identity for this operation, here’s how to get them

without triggering a warning message to stderr.

>>> x[:,1] = quiet_log(0.0)
>>> x
array([[-3.2, -inf],
       [-3.2, -inf],
       [-3.2, -inf]])

>>> logsumexp_array(x, axis=1)
array([-3.2, -3.2, -3.2])

>>> w = np.zeros((3,3), dtype=float) - 3.2
>>> logsumexp(w, w)
array([[-2.50685282, -2.50685282, -2.50685282],
       [-2.50685282, -2.50685282, -2.50685282],
       [-2.50685282, -2.50685282, -2.50685282]])

>>> z = log_zeros((3,3), dtype=float)
>>> (logsumexp(w, z) == w).all()
True

>>> y = np.zeros((2,3), dtype=float) - 20.0
>>> y[:,1] = quiet_log(0.0)
>>> y
array([[-20., -inf, -20.],
       [-20., -inf, -20.]])

>>> log_dot(x,y)
array([[-23.2,  -inf, -23.2],
       [-23.2,  -inf, -23.2],
       [-23.2,  -inf, -23.2]])

>>> log_dot(y,x)
array([[-22.50685282,         -inf],
       [-22.50685282,         -inf]])

>>> x_exp = np.exp(x)
>>> y_exp = np.exp(y)
>>> dot_exp = np.dot(y_exp, x_exp)
>>> (log_dot(y,x) == quiet_log(dot_exp)).all()
True

Note that using a single maximum value can run into underflow problems pretty easily if the values in the array are very different from each other (about 705 for ordinary doubles). Using an array of bias values makes this work more correctly.

>>> test_arr = np.array((-1.0, -7006.1, -7260.1, -7469.1, -7568.2, -7675.8, -7814.9, -9617.2), dtype = float)
>>> test_arr[0] = quiet_log(0.0)
>>> print test_arr
[   -inf -7006.1 -7260.1 -7469.1 -7568.2 -7675.8 -7814.9 -9617.2]
>>> z = log_zeros((len(test_arr),), dtype=float)
>>> y = logsumexp(test_arr, z)
>>> print y
[   -inf -7006.1 -7260.1 -7469.1 -7568.2 -7675.8 -7814.9 -9617.2]
>>> (y == test_arr).all()
True

>>> test_arr2 = test_arr - 10000.0
>>> test_arr2
array([    -inf, -17006.1, -17260.1, -17469.1, -17568.2, -17675.8,
       -17814.9, -19617.2])

As might be expected, ‘adding’ numbers which are so much closer to ‘0’ makes no difference. In fact, the values added were actually 0 in this example.

>>> logsumexp(test_arr, test_arr2)
array([   -inf, -7006.1, -7260.1, -7469.1, -7568.2, -7675.8, -7814.9,
       -9617.2])

This module includes some handy constants:

>>> floatutils.float_to_readable_string(RECIP_LOG2)
'+(+0000)0x71547652b82fe'

>>> floatutils.float_to_readable_string(LOG2_TO_LOG10_FACTOR)
'+(-0002)0x34413509f79fe'

class onyx.util.mathutils.Primes

Bases: object

Generates primes based on the Eratosthenes sieve, but not requiring an end limit. Instance is intended to be used (and reused) as an iterator or an indexed container.

Module attribute ‘primes’ is an instance of Primes() intended for general use.

>>> primes[20]
73

>>> tuple(islice(primes, 100, 101))
(547,)
>>> primes[42:44]
[191, 193]

>>> primes[:]
Traceback (most recent call last):
   ...
ValueError: expected non-negative value for indexing or slice.stop, got None

extend()

onyx.util.mathutils.distance(v1, v2)

Euclidean distance between two vectors

>>> orig = np.zeros(4)
>>> one = np.ones(4)
>>> d = distance(orig, one)
>>> floatutils.float_to_readable_string(d)
'+(+0001)0x0000000000000'
>>> floatutils.float_to_readable_string(distance(orig, 2 * one))
'+(+0002)0x0000000000000'

onyx.util.mathutils.factors(value)

Return the ordered sequence of the factors of positive integer ‘value’.

See also prime_factors().

>>> factors(720)
(1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 30, 36, 40, 45, 48, 60, 72, 80, 90, 120, 144, 180, 240, 360, 720)

>>> for i in xrange(1, 13):
...   i, factors(i)
(1, (1,))
(2, (1, 2))
(3, (1, 3))
(4, (1, 2, 4))
(5, (1, 5))
(6, (1, 2, 3, 6))
(7, (1, 7))
(8, (1, 2, 4, 8))
(9, (1, 3, 9))
(10, (1, 2, 5, 10))
(11, (1, 11))
(12, (1, 2, 3, 4, 6, 12))

>>> big = 2 * 2 * 3 * 3 * 5 * 5 * 7 * 7 * 11 * 11 * 13 * 13
>>> big
901800900
>>> f = factors(big)
>>> len(f)
729

One invariant of the factors

>>> tuple((f1, f2, f1*f2) for f1, f2 in izip(f, reversed(f)) if f1 * f2 != big)
()

onyx.util.mathutils.find_cumulative_index(arr, limit)

Find the index of the first element in an iterable of values whose cumulative value is at or above a given limit

Returns the index, or None if the total of the entire vector is greater than the limit.

>>> arr = np.array((1, 2, 3, 4))
>>> find_cumulative_index(arr, 0)
0
>>> find_cumulative_index(arr, 3)
1
>>> find_cumulative_index(arr, 7)
3
>>> find_cumulative_index(arr, 11) is None
True

onyx.util.mathutils.log_dot(a, b)

Compute the dot product of two 2-dimensional Numpy arrays in the log domain.

onyx.util.mathutils.log_zeros(shape, dtype=<type 'float'>)

Produce an array with the given shape filled with log-domain 0 values, i.e., -inf.

onyx.util.mathutils.logsumexp(x, y)

Compute log(exp(x) + exp(y))) where x and y are either scalars, Numpy arrays of the same shape, or Numpy arrays that can be projected into the same shape.

onyx.util.mathutils.logsumexp_array(x, axis=None)

Compute log(sum(exp(x))) along a particular axis for Numpy arrays.

onyx.util.mathutils.main(args)

onyx.util.mathutils.prime_factors(value)

Return the ordered sequence of the prime factors of non-negative integer ‘value’.

See also factors().

Look at some factorizations

>>> for i in xrange(20):
...   i, prime_factors(i)
(0, [0])
(1, [1])
(2, [2])
(3, [3])
(4, [2, 2])
(5, [5])
(6, [2, 3])
(7, [7])
(8, [2, 2, 2])
(9, [3, 3])
(10, [2, 5])
(11, [11])
(12, [2, 2, 3])
(13, [13])
(14, [2, 7])
(15, [3, 5])
(16, [2, 2, 2, 2])
(17, [17])
(18, [2, 3, 3])
(19, [19])

Verify, for a few factorizations, that the product of the returned items equals the argument

>>> tuple((i, prime_factors(i)) for i in xrange(0, 1000) if reduce(mul, prime_factors(i)) != i)
()

Errors

>>> prime_factors(-1)
Traceback (most recent call last):
...
ValueError: expected a non-negative number, got -1

>>> prime_factors('a')
Traceback (most recent call last):
...
TypeError: expected positive 'int', got 'str'

onyx.util.mathutils.quiet_log(item)

Return a copy of item that’s been transformed to the log-domain, but do the operation quietly even if there are zeros in the item.

onyx.util.mathutils.rollaxis(a, axis, start=0)

Roll the specified axis until it lies in a given position.

Parameters

a : ndarray

Input array.

axis : int

The axis to roll. The positions of the other axes do not change relative to one another.

start : int, optional

The axis is rolled until it lies before this position.

Returns

res : ndarray

Output array.

Examples

>>> a = np.ones((3,4,5,6))
>>> rollaxis(a, 3, 1).shape
(3, 6, 4, 5)
>>> rollaxis(a, 2).shape
(5, 3, 4, 6)
>>> rollaxis(a, 1, 4).shape
(3, 5, 6, 4)

First axis becomes the last:

>>> rollaxis(a, 0, -1).shape
(4, 5, 6, 3)

Last axis becomes the first:

>>> rollaxis(a, -1).shape
(6, 3, 4, 5)

onyx.util.mathutils.safe_log_divide(x, y)

Compute x - y where x and y are either scalars, Numpy arrays of the same shape, or Numpy arrays that can be projected into the same shape, and make sure that (-inf) - (-inf) is handled correctly and quietly. Complain if the second operand is -inf but the first is not.

>>> zero = quiet_log(0.0)
>>> result = safe_log_divide(zero, 1.0)
>>> result == zero
True

>>> safe_log_divide(1.0, zero)
Traceback (most recent call last):
ValueError: log division by zero

>>> result = safe_log_divide(zero, zero)
>>> result == zero
True

>>> x = np.zeros((3,2), dtype=float) - 3.2
>>> x[:,1] = quiet_log(0.0)
>>> x
array([[-3.2, -inf],
       [-3.2, -inf],
       [-3.2, -inf]])

>>> safe_log_divide(x, x)
array([[  0., -inf],
       [  0., -inf],
       [  0., -inf]])

>>> y = np.zeros((3,2), dtype=float) - 3.2
>>> safe_log_divide(x, y)
array([[  0., -inf],
       [  0., -inf],
       [  0., -inf]])

>>> safe_log_divide(y, x)
Traceback (most recent call last):
ValueError: log division by zero

onyx.util.mathutils.truncate_fp_array_in_place(array)

>>> arr0 = np.arange(16, dtype=np.float64)
>>> arr0 += 1 + (1<<6) + (1<<28)
>>> arr0 = 1 / arr0
>>> for x in arr0: print floatutils.float_to_readable_string(x)
+(-0029)0xfffff7e000210
+(-0029)0xfffff7c000220
+(-0029)0xfffff7a000231
+(-0029)0xfffff78000242
+(-0029)0xfffff76000253
+(-0029)0xfffff74000264
+(-0029)0xfffff72000276
+(-0029)0xfffff70000288
+(-0029)0xfffff6e00029a
+(-0029)0xfffff6c0002ac
+(-0029)0xfffff6a0002bf
+(-0029)0xfffff680002d2
+(-0029)0xfffff660002e5
+(-0029)0xfffff640002f8
+(-0029)0xfffff6200030c
+(-0029)0xfffff60000320

>>> truncate_fp_array_in_place(arr0)
>>> for x in arr0: print floatutils.float_to_readable_string(x)
+(-0029)0xfffff70000000
+(-0029)0xfffff70000000
+(-0029)0xfffff70000000
+(-0029)0xfffff70000000
+(-0029)0xfffff70000000
+(-0029)0xfffff70000000
+(-0029)0xfffff70000000
+(-0029)0xfffff70000000
+(-0029)0xfffff60000000
+(-0029)0xfffff60000000
+(-0029)0xfffff60000000
+(-0029)0xfffff60000000
+(-0029)0xfffff60000000
+(-0029)0xfffff60000000
+(-0029)0xfffff60000000
+(-0029)0xfffff60000000