From 069122cd7eb8cee87903a280aac2ca57456efd85 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Mon, 29 Feb 2016 09:34:56 +0100 Subject: [PATCH 01/42] Moved util files to more appropriate package --- accumarray.py | 108 +++++++++++++++++++++++++++++++++++++++++++++++ setmembership.py | 25 +++++++++++ 2 files changed, 133 insertions(+) create mode 100644 accumarray.py create mode 100644 setmembership.py diff --git a/accumarray.py b/accumarray.py new file mode 100644 index 0000000000..84bc1caa88 --- /dev/null +++ b/accumarray.py @@ -0,0 +1,108 @@ +from itertools import product +import numpy as np + + +def accum(accmap, a, func=None, size=None, fill_value=0, dtype=None): + """ + An accumulation function similar to Matlab's `accumarray` function. + + Parameters + ---------- + accmap : ndarray + This is the "accumulation map". It maps input (i.e. indices into + `a`) to their destination in the output array. The first `a.ndim` + dimensions of `accmap` must be the same as `a.shape`. That is, + `accmap.shape[:a.ndim]` must equal `a.shape`. For example, if `a` + has shape (15,4), then `accmap.shape[:2]` must equal (15,4). In this + case `accmap[i,j]` gives the index into the output array where + element (i,j) of `a` is to be accumulated. If the output is, say, + a 2D, then `accmap` must have shape (15,4,2). The value in the + last dimension give indices into the output array. If the output is + 1D, then the shape of `accmap` can be either (15,4) or (15,4,1) + a : ndarray + The input data to be accumulated. + func : callable or None + The accumulation function. The function will be passed a list + of values from `a` to be accumulated. + If None, numpy.sum is assumed. + size : ndarray or None + The size of the output array. If None, the size will be determined + from `accmap`. + fill_value : scalar + The default value for elements of the output array. + dtype : numpy data type, or None + The data type of the output array. If None, the data type of + `a` is used. + + Returns + ------- + out : ndarray + The accumulated results. + + The shape of `out` is `size` if `size` is given. Otherwise the + shape is determined by the (lexicographically) largest indices of + the output found in `accmap`. + + + Examples + -------- + >>> from numpy import array, prod + >>> a = array([[1,2,3],[4,-1,6],[-1,8,9]]) + >>> a + array([[ 1, 2, 3], + [ 4, -1, 6], + [-1, 8, 9]]) + >>> # Sum the diagonals. + >>> accmap = array([[0,1,2],[2,0,1],[1,2,0]]) + >>> s = accum(accmap, a) + array([9, 7, 15]) + >>> # A 2D output, from sub-arrays with shapes and positions like this: + >>> # [ (2,2) (2,1)] + >>> # [ (1,2) (1,1)] + >>> accmap = array([ + [[0,0],[0,0],[0,1]], + [[0,0],[0,0],[0,1]], + [[1,0],[1,0],[1,1]], + ]) + >>> # Accumulate using a product. + >>> accum(accmap, a, func=prod, dtype=float) + array([[ -8., 18.], + [ -8., 9.]]) + >>> # Same accmap, but create an array of lists of values. + >>> accum(accmap, a, func=lambda x: x, dtype='O') + array([[[1, 2, 4, -1], [3, 6]], + [[-1, 8], [9]]], dtype=object) + """ + + # Check for bad arguments and handle the defaults. + if accmap.shape[:a.ndim] != a.shape: + raise ValueError("The initial dimensions of accmap must be the same as a.shape") + if func is None: + func = np.sum + if dtype is None: + dtype = a.dtype + if accmap.shape == a.shape: + accmap = np.expand_dims(accmap, -1) + adims = tuple(range(a.ndim)) + if size is None: + size = 1 + np.squeeze(np.apply_over_axes(np.max, accmap, axes=adims)) + size = np.atleast_1d(size) + + # Create an array of python lists of values. + vals = np.empty(size, dtype='O') + for s in product(*[range(k) for k in size]): + vals[s] = [] + for s in product(*[range(k) for k in a.shape]): + indx = tuple(accmap[s]) + val = a[s] + vals[indx].append(val) + + # Create the output array. + out = np.empty(size, dtype=dtype) + for s in product(*[range(k) for k in size]): + if vals[s] == []: + out[s] = fill_value + else: + out[s] = func(vals[s]) + + return out diff --git a/setmembership.py b/setmembership.py new file mode 100644 index 0000000000..1467f76ada --- /dev/null +++ b/setmembership.py @@ -0,0 +1,25 @@ +# -*- coding: utf-8 -*- +""" +Created on Thu Feb 25 20:16:42 2016 + +@author: keile +""" + +import numpy as np + + +def unique_rows(data): + """ + Function similar to Matlab's unique(...,'rows') + + Copied from + http://stackoverflow.com/questions/16970982/find-unique-rows-in-numpy-array/ + (summary pretty far down on the page) + Note: I have no idea what happens here + """ + b = np.ascontiguousarray(data).view(np.dtype((np.void, + data.dtype.itemsize * data.shape[1]))) + _, ia = np.unique(b, return_index=True) + _, ic = np.unique(b, return_inverse=True) +# return np.unique(b).view(data.dtype).reshape(-1, data.shape[1]), ia, ic + return data[ia], ia, ic From 6cea356960a41946dd709e009f63d88e729454a3 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Tue, 1 Mar 2016 20:05:01 +0100 Subject: [PATCH 02/42] Methods to compress and recover matrices (rlencode and -decode) --- matrix_compression.py | 48 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 48 insertions(+) create mode 100644 matrix_compression.py diff --git a/matrix_compression.py b/matrix_compression.py new file mode 100644 index 0000000000..12e6efff81 --- /dev/null +++ b/matrix_compression.py @@ -0,0 +1,48 @@ +# -*- coding: utf-8 -*- +""" +Created on Tue Mar 1 08:36:47 2016 + +@author: keile +""" + +import numpy as np + +def rldecode(A,n): + """ Decode compressed information. + + The code is heavily inspired by MRST's function with the same name, + however, requirements on the shape of functions are probably somewhat + different. + + >>> rlencode(np.array([1, 2, 3]), np.array([2, 3, 1])) + [1, 1, 2, 2, 2, 3] + + >>> rlencode(np.array([0, 2]), np.array([0, 3])) + [2, 2, 2] + + Args: + A (double, m x k), compressed matrix to be recovered. The + compression should be along dimension 1 + n (int): Number of occurences for each element + """ + r = n > 0 + i = np.cumsum(np.hstack((np.zeros(1), n[r])), dtype='>i4') + j = np.zeros(i[-1]) + j[i[1:-1:]] = 1 + B = A[np.cumsum(j, dtype='>i4')] + return B + +def rlencode(A): + """ Compress matrix by looking for identical columns. """ + comp = A[::, 0:-1] != A[::, 1::] + i = np.any(comp, axis=0) + i = np.hstack((np.argwhere(i).ravel(), (A.shape[1]-1))) + + num = np.diff(np.hstack((np.array([-1]), i))) + + return A[::, i], num + + +if __name__ == '__main__': + A = np.array([[1, 2, 2, 3],[2, 2, 2, 3]]) + rlencode(A) \ No newline at end of file From d71232a0f49160b0c4c006ad04d3892e0fe36247 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Tue, 1 Mar 2016 22:00:40 +0100 Subject: [PATCH 03/42] Bugfix in doctests in matrix_compression (rlencode decode) --- matrix_compression.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/matrix_compression.py b/matrix_compression.py index 12e6efff81..1477fa3862 100644 --- a/matrix_compression.py +++ b/matrix_compression.py @@ -14,11 +14,11 @@ def rldecode(A,n): however, requirements on the shape of functions are probably somewhat different. - >>> rlencode(np.array([1, 2, 3]), np.array([2, 3, 1])) + >>> rldecode(np.array([1, 2, 3]), np.array([2, 3, 1])) [1, 1, 2, 2, 2, 3] - >>> rlencode(np.array([0, 2]), np.array([0, 3])) - [2, 2, 2] + >>> rldecode(np.array([1, 2]), np.array([1, 3])) + [1, 2, 2, 2] Args: A (double, m x k), compressed matrix to be recovered. The From e368e3b07d07ff0d263c67374725bba104026dce Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Tue, 19 Apr 2016 09:39:13 +0200 Subject: [PATCH 04/42] Emulate behavior of MRST mcolon --- mcolon.py | 43 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 43 insertions(+) create mode 100644 mcolon.py diff --git a/mcolon.py b/mcolon.py new file mode 100644 index 0000000000..6309e984d1 --- /dev/null +++ b/mcolon.py @@ -0,0 +1,43 @@ +# -*- coding: utf-8 -*- +""" +Created on Fri Mar 4 21:02:56 2016 + +@author: keile +""" + +import numpy as np + + +def mcolon(lo, hi): + """ Get set of expanded indices + + >>> mcolon(np.array([0, 0, 0]), np.array([1, 3, 2])) + array([0, 1, 0, 1, 2, 3, 0, 1, 2]) + + >>> mcolon(np.array([0, 1]), np.array([2])) + array([0, 1, 2, 1, 2]) + + >>> mcolon(np.array([0, 1, 1, 1]), np.array([0, 3, 3, 3])) + array([0, 1, 2, 3, 1, 2, 3, 1, 2, 3]) + """ + if lo.size == 1: + lo = lo * np.ones(hi.size, dtype='int64') + if hi.size == 1: + hi = hi * np.ones(lo.size, dtype='int64') + if lo.size != hi.size: + raise ValueError('Low and high should have same number of elements, ' + 'or a single item ') + + i = hi >= lo + if not any(i): + return None + + lo = lo[i] + hi = hi[i] + d = hi - lo + 1 + n = np.sum(d) + + x = np.ones(n, dtype='int64') + x[0] = lo[0] + x[np.cumsum(d[0:-1]).astype('int64')] = lo[1:] - hi[0:-1] + return np.cumsum(x) From 85d89eaf910457b8a42be0352829a4b31a572e90 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Wed, 4 May 2016 12:16:36 +0200 Subject: [PATCH 05/42] Upload __init__.py files. Should have done this a long time ago --- __init__.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 __init__.py diff --git a/__init__.py b/__init__.py new file mode 100644 index 0000000000..e69de29bb2 From 9f64225e7f19eb5ff410c89fcd501f79eacdbe80 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Tue, 31 May 2016 11:12:07 +0200 Subject: [PATCH 06/42] Added equivalent of matlab ismember_rows. Implementation not yet streamlined. --- setmembership.py | 79 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 79 insertions(+) diff --git a/setmembership.py b/setmembership.py index 1467f76ada..ac15068a32 100644 --- a/setmembership.py +++ b/setmembership.py @@ -23,3 +23,82 @@ def unique_rows(data): _, ic = np.unique(b, return_inverse=True) # return np.unique(b).view(data.dtype).reshape(-1, data.shape[1]), ia, ic return data[ia], ia, ic + + +def asvoid(arr): + """ + + Taken from + http://stackoverflow.com/questions/22699756/python-version-of-ismember-with-rows-and-index + + View the array as dtype np.void (bytes) + This views the last axis of ND-arrays as bytes so you can perform + comparisons on the entire row. + http://stackoverflow.com/a/16840350/190597 (Jaime, 2013-05) + Warning: When using asvoid for comparison, note that float zeros may + compare UNEQUALLY + >>> asvoid([-0.]) == asvoid([0.]) + array([False], dtype=bool) + """ + arr = np.ascontiguousarray(arr) + return arr.view(np.dtype((np.void, arr.dtype.itemsize * arr.shape[-1]))) + + +def find_occ(a, b): + """ + + + Code snippet found at + http://stackoverflow.com/questions/15864082/python-equivalent-of-matlabs-ismember-function?rq=1 + + Parameters + ---------- + a + b + + Returns + ------- + + """ + bind = {} + for i, elt in enumerate(b): + if elt not in bind: + bind[elt] = i + return [bind.get(itm, None) for itm in a] + + +def ismember_rows(a, b, sort=True): + """ + Examples + >>> a = np.array([[1, 3, 3, 1, 7], [3, 3, 2, 3, 0]]) + >>> b = np.array([[3, 1, 3, 5, 3], [3, 3, 2, 1, 2]]) + >>> ismember_rows(a, b) + (array([ True, True, True, True, False], dtype=bool), [1, 0, 2, 1]) + + >>> a = np.array([[1, 3, 3, 1, 7], [3, 3, 2, 3, 0]]) + >>> b = np.array([[3, 1, 2, 5, 1], [3, 3, 3, 1, 2]]) + >>> ismember_rows(a, b, sort=False) + (array([ True, True, False, True, False], dtype=bool), [1, 0, 1]) + """ + if sort: + sa = np.sort(a, axis=0) + sb = np.sort(b, axis=0) + else: + sa = a + sb = b + + voida = asvoid(sa.transpose()) + voidb = asvoid(sb.transpose()) + unq, j, k, count = np.unique(np.vstack((voida, voidb)), return_index=True, + return_inverse=True, return_counts=True) + + num_a = a.shape[1] + ind_a = np.arange(num_a) + ind_b = num_a + np.arange(b.shape[1]) + num_occ_a = count[k[ind_a]] + ismem_a = (num_occ_a > 1) + occ_a = k[ind_a[ismem_a]] + occ_b = k[ind_b] + + ind_of_a_in_b = find_occ(occ_a, occ_b) + return ismem_a, ind_of_a_in_b From b4653915ae84c99143927b61fa9fd6c73a283b49 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Tue, 7 Jun 2016 10:19:53 +0200 Subject: [PATCH 07/42] Fix of doctest in utils.mcolon --- mcolon.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/mcolon.py b/mcolon.py index 6309e984d1..1eef9cb949 100644 --- a/mcolon.py +++ b/mcolon.py @@ -12,13 +12,13 @@ def mcolon(lo, hi): """ Get set of expanded indices >>> mcolon(np.array([0, 0, 0]), np.array([1, 3, 2])) - array([0, 1, 0, 1, 2, 3, 0, 1, 2]) + array([0, 1, 0, 1, 2, 3, 0, 1, 2], dtype=int64) >>> mcolon(np.array([0, 1]), np.array([2])) - array([0, 1, 2, 1, 2]) + array([0, 1, 2, 1, 2], dtype=int64) >>> mcolon(np.array([0, 1, 1, 1]), np.array([0, 3, 3, 3])) - array([0, 1, 2, 3, 1, 2, 3, 1, 2, 3]) + array([0, 1, 2, 3, 1, 2, 3, 1, 2, 3], dtype=int64) """ if lo.size == 1: lo = lo * np.ones(hi.size, dtype='int64') From 962eaea3b7a9fb3c911e2493486aa3de269387dd Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Fri, 21 Oct 2016 13:16:07 +0200 Subject: [PATCH 08/42] Move unit tests for utility functions --- tests/utils/__init__.py | 0 tests/utils/test_mcolon.py | 24 +++++++++++++++ tests/utils/test_setmembership.py | 50 +++++++++++++++++++++++++++++++ 3 files changed, 74 insertions(+) create mode 100644 tests/utils/__init__.py create mode 100644 tests/utils/test_mcolon.py create mode 100644 tests/utils/test_setmembership.py diff --git a/tests/utils/__init__.py b/tests/utils/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/tests/utils/test_mcolon.py b/tests/utils/test_mcolon.py new file mode 100644 index 0000000000..dade528d2d --- /dev/null +++ b/tests/utils/test_mcolon.py @@ -0,0 +1,24 @@ +# -*- coding: utf-8 -*- +""" +Created on Thu Apr 14 13:51:05 2016 + +@author: eke001 +""" + +import numpy as np + +from utils import mcolon + + +def test_mcolon_simple(): + a = np.array([1, 2]) + b = np.array([2, 3]) + c = mcolon.mcolon(a, b) + assert np.all((c - np.array([1, 2, 2, 3])) == 0) + + +def test_mcolon_one_missing(): + a = np.array([1, 2]) + b = np.array([2, 1]) + c = mcolon.mcolon(a, b) + assert np.all((c - np.array([1, 2])) == 0) diff --git a/tests/utils/test_setmembership.py b/tests/utils/test_setmembership.py new file mode 100644 index 0000000000..94a9d11709 --- /dev/null +++ b/tests/utils/test_setmembership.py @@ -0,0 +1,50 @@ +# -*- coding: utf-8 -*- +""" +Created on Thu Feb 25 20:43:38 2016 + +@author: keile +""" +import numpy as np + +from utils import setmembership + +def test_unique_rows_1(): + + a = np.array([[1, 2], [2, 1], [2, 4], [2, 1], [2, 4]]) + ua_expected = np.array([[1, 2], [2, 1], [2, 4]]) + ia_expected = np.array([0, 1, 2]) + ic_expected = np.array([0, 1, 2, 1, 2]) + ua, ia, ic = setmembership.unique_rows(a) + assert np.sum(np.abs(ua) - np.abs(ua_expected)) == 0 + assert np.all(ia - ia_expected == 0) + assert np.all(ic - ic_expected == 0) + + +def test_ismember_rows_with_sort(): + a = np.array([[1, 3, 3, 1, 7], [3, 3, 2, 3, 0]]) + b = np.array([[3, 1, 3, 5, 3], [3, 3, 2, 1, 2]]) + ma, ia = setmembership.ismember_rows(a, b) + + ma_known = np.array([1, 1, 1, 1, 0], dtype=bool) + ia_known = np.array([1, 0, 2, 1]) + + assert np.allclose(ma, ma_known) + assert np.allclose(ia, ia_known) + + +def test_ismember_rows_no_sort(): + a = np.array([[1, 3, 3, 1, 7], [3, 3, 2, 3, 0]]) + b = np.array([[3, 1, 2, 5, 3], [3, 3, 3, 1, 1]]) + ma, ia = setmembership.ismember_rows(a, b, sort=False) + + ma_known = np.array([1, 1, 0, 1, 0], dtype=bool) + ia_known = np.array([1, 0, 1]) + + assert np.allclose(ma, ma_known) + assert np.allclose(ia, ia_known) + + +if __name__ == '__main__': + test_unique_rows_1() + test_ismember_rows_with_sort() + test_ismember_rows_no_sort() \ No newline at end of file From a15695e299fbb773c54a9ba5c0b40e58bc3805bd Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Mon, 19 Dec 2016 09:30:01 +0100 Subject: [PATCH 09/42] Gitignore --- .gitignore | 1 + 1 file changed, 1 insertion(+) create mode 100644 .gitignore diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000000..0d20b6487c --- /dev/null +++ b/.gitignore @@ -0,0 +1 @@ +*.pyc From c5ed565b2ad1067500dbac326e0b4f1b00d9a488 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Mon, 2 Jan 2017 20:14:21 +0100 Subject: [PATCH 10/42] Moved unit tests to more appropriate folder. --- tests/__init__.py | 0 tests/{utils => }/test_mcolon.py | 0 tests/{utils => }/test_setmembership.py | 0 3 files changed, 0 insertions(+), 0 deletions(-) create mode 100644 tests/__init__.py rename tests/{utils => }/test_mcolon.py (100%) rename tests/{utils => }/test_setmembership.py (100%) diff --git a/tests/__init__.py b/tests/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/tests/utils/test_mcolon.py b/tests/test_mcolon.py similarity index 100% rename from tests/utils/test_mcolon.py rename to tests/test_mcolon.py diff --git a/tests/utils/test_setmembership.py b/tests/test_setmembership.py similarity index 100% rename from tests/utils/test_setmembership.py rename to tests/test_setmembership.py From 4dffb6ac1b13a2ba40255e0910600da8a57d5379 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Thu, 5 Jan 2017 16:27:59 +0100 Subject: [PATCH 11/42] Removed __init__ from unused directory --- tests/utils/__init__.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) delete mode 100644 tests/utils/__init__.py diff --git a/tests/utils/__init__.py b/tests/utils/__init__.py deleted file mode 100644 index e69de29bb2..0000000000 From 905fb69d33261ffdbe05cac3f1924919ed33221f Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Mon, 9 Jan 2017 11:31:48 +0100 Subject: [PATCH 12/42] Generator for multinary numbers of given base and length. Also unit tests. --- permutations.py | 46 ++++++++++++++++++++++++++++++++++++++ tests/test_permutations.py | 46 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 92 insertions(+) create mode 100644 permutations.py create mode 100644 tests/test_permutations.py diff --git a/permutations.py b/permutations.py new file mode 100644 index 0000000000..48aec5b164 --- /dev/null +++ b/permutations.py @@ -0,0 +1,46 @@ + +def multinary_permutations(base, length): + """ + Define a generator over all numbers of a certain length for a number system + with a specified base. + + For details on the decomposition into an arbitrary base see + http://math.stackexchange.com/questions/111150/changing-a-number-between-arbitrary-bases + + Note that the generator will loop over base**length combinations. + + Examples: + + Construct the numbers [00] to [11] in binary numbers + >>> multinary_permutations(2, 2) + [array([ 0., 0.]), array([ 1., 0.]), array([ 0., 1.]), array([ 1., 1.])] + + Construct the numbers from 0 to 99 (permuted) in the decimal number + system. + >>> it = multinary_permutation(10, 2) + + Parameters: + base (int): Base of the number system + length (int): Number of digits in the numbers + + Yields: + array, size length: Array describing the next number combination. + + """ + + # There are in total base ** length numbers to be covered, these need to be + # rewritten into the base number system + for iter1 in range(base ** length): + + # Array to store the multi-d index of the current index + bit_val = [0] * length + # Number to be decomposed + v = iter1 + + # Loop over all digits, find the expression of v in that system + for iter2 in range(length): + bit_val[iter2] = v % base + v = v // base + # Yield the next value + yield bit_val + diff --git a/tests/test_permutations.py b/tests/test_permutations.py new file mode 100644 index 0000000000..67359d5f46 --- /dev/null +++ b/tests/test_permutations.py @@ -0,0 +1,46 @@ +import unittest +import numpy as np + +from utils import permutations + +class TestPermutations(unittest.TestCase): + + def compare_lists(self, base, length, lst): + # Compare a pre-defined list with the result of multinary_permutations + # Define a generator, and check that all values produced are contained within lst + # Also count the number of iterations + iter_cnt = 0 + for a in permutations.multinary_permutations(base, length): + found = False + for b in lst: + if np.array_equal(np.array(a), np.array(b)): + found = True + break + assert found + iter_cnt += 1 + assert iter_cnt == len(lst) + + + def test_length_2(self): + # Explicitly construct a 2D array of all combination of base 3 + base = 3 + length = 2 + lst = [] + for i in range(base): + for j in range(base): + lst.append([i, j]) + self.compare_lists(base, length, lst) + + def test_base_4(self): + # Check with a manually constructed list of length 3 + base = 4 + length = 3 + lst = [] + for i in range(base): + for j in range(base): + for k in range(base): + lst.append([i, j, k]) + self.compare_lists(base, length, lst) + + if __name__ == '__main__': + unittest.main() From adee794550a258c9307eb0d9e343840abef1bd2f Mon Sep 17 00:00:00 2001 From: Runar Date: Thu, 2 Feb 2017 14:30:59 +0100 Subject: [PATCH 13/42] New function: half_space_int - find points that lie in the intersection of half planes Also with unit tests --- .gitignore | 2 ++ half_space_int.py | 48 ++++++++++++++++++++++++++++++++++++ tests/test_half_space_int.py | 26 +++++++++++++++++++ 3 files changed, 76 insertions(+) create mode 100644 half_space_int.py create mode 100644 tests/test_half_space_int.py diff --git a/.gitignore b/.gitignore index 0d20b6487c..840e386596 100644 --- a/.gitignore +++ b/.gitignore @@ -1 +1,3 @@ *.pyc +*.py# +*.py~ \ No newline at end of file diff --git a/half_space_int.py b/half_space_int.py new file mode 100644 index 0000000000..23f3f0d15c --- /dev/null +++ b/half_space_int.py @@ -0,0 +1,48 @@ +import numpy as np + + +def half_space_int(n,x0, pts): + """ + Find the points that lie in the intersection of half spaces (3D) + + Parameters + ---------- + n : ndarray + This is the normal vectors of the half planes. The normal + vectors is assumed to point out of the half spaces. + x0 : ndarray + Point on the boundary of the half-spaces. Half space i is given + by all points satisfying (x - x0[:,i])*n[:,i]<=0 + pts : ndarray + The points to be tested if they are in the intersection of all + half-spaces or not. + + Returns + ------- + out : ndarray + A logical array with length equal number of pts. + + out[i] is True if pts[:,i] is in all half-spaces + + + Examples + -------- + >>> import numpy as np + >>> n = np.array([[0,1],[1,0],[0,0]]) + >>> x0 = np.array([[0,-1],[0,0],[0,0]]) + >>> pts = np.array([[-1,1,4],[2,1,-2],[0,0,0]]) + >>> out = half_plane_int(n,x0,pts) + array([False, True, False], dtype=bool) + """ + assert n.shape[0] == 3, ' only 3D supported' + assert x0.shape[0] == 3, ' only 3D supported' + assert pts.shape[0] == 3, ' only 3D supported' + assert n.shape[1] == x0.shape[1], 'ther must be the same number of normal vectors as points' + n_pts = pts.shape[1] + in_hull = np.zeros(n_pts) + x0 = np.repeat(x0[:,:,np.newaxis], n_pts,axis=2) + n = np.repeat( n[:,:,np.newaxis], n_pts,axis=2) + for i in range(x0.shape[1]): + in_hull += np.sum((pts - x0[:,i,:])*n[:,i,:],axis=0)<=0 + + return in_hull==x0.shape[1] diff --git a/tests/test_half_space_int.py b/tests/test_half_space_int.py new file mode 100644 index 0000000000..8887917763 --- /dev/null +++ b/tests/test_half_space_int.py @@ -0,0 +1,26 @@ +import unittest +import numpy as np + +from utils import half_space_int + + +def test_one_half_space(): + + n = np.array([[-1],[0],[0]]) + x0 = np.array([[0],[0],[0]]) + pts = np.array([[1,-1],[0,0],[0,0]]) + out = half_space_int.half_space_int(n,x0,pts) + assert np.all(out == np.array([True,False])) + +def test_two_half_spaces(): + n = np.array([[-1,0],[0,-1],[0,0]]) + x0 = np.array([[0,0],[0,1],[0,0]]) + pts = np.array([[1,-1,1,0],[2,0,2,0],[0,0,0,0]]) + out = half_space_int.half_space_int(n,x0,pts) + assert np.all(out == np.array([True,False,True,False])) + + +if __name__ == '__main__': + test_one_half_space() + test_two_half_spaces() + From 23dc8c34229e8d3deefdad95daa798fb3d99d3cb Mon Sep 17 00:00:00 2001 From: Runar Date: Thu, 2 Feb 2017 16:13:07 +0100 Subject: [PATCH 14/42] Fixed a mistake in the example in half_space_int --- half_space_int.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/half_space_int.py b/half_space_int.py index 23f3f0d15c..f5ad0f7e59 100644 --- a/half_space_int.py +++ b/half_space_int.py @@ -30,8 +30,8 @@ def half_space_int(n,x0, pts): >>> import numpy as np >>> n = np.array([[0,1],[1,0],[0,0]]) >>> x0 = np.array([[0,-1],[0,0],[0,0]]) - >>> pts = np.array([[-1,1,4],[2,1,-2],[0,0,0]]) - >>> out = half_plane_int(n,x0,pts) + >>> pts = np.array([[-1,-1,4],[2,-2,-2],[0,0,0]]) + >>> half_space_int(n,x0,pts) array([False, True, False], dtype=bool) """ assert n.shape[0] == 3, ' only 3D supported' From c4ddc6d6bc983cf2db546cdd100b64314863a9f7 Mon Sep 17 00:00:00 2001 From: Alessio Fumagalli Date: Fri, 3 Feb 2017 16:07:02 +0100 Subject: [PATCH 15/42] use the standard unit test procedure --- tests/test_half_space_int.py | 33 ++++++++++++++++++--------------- 1 file changed, 18 insertions(+), 15 deletions(-) diff --git a/tests/test_half_space_int.py b/tests/test_half_space_int.py index 8887917763..c4a7342663 100644 --- a/tests/test_half_space_int.py +++ b/tests/test_half_space_int.py @@ -3,24 +3,27 @@ from utils import half_space_int +#------------------------------------------------------------------------------# -def test_one_half_space(): +class BasicsTest( unittest.TestCase ): - n = np.array([[-1],[0],[0]]) - x0 = np.array([[0],[0],[0]]) - pts = np.array([[1,-1],[0,0],[0,0]]) - out = half_space_int.half_space_int(n,x0,pts) - assert np.all(out == np.array([True,False])) +#------------------------------------------------------------------------------# -def test_two_half_spaces(): - n = np.array([[-1,0],[0,-1],[0,0]]) - x0 = np.array([[0,0],[0,1],[0,0]]) - pts = np.array([[1,-1,1,0],[2,0,2,0],[0,0,0,0]]) - out = half_space_int.half_space_int(n,x0,pts) - assert np.all(out == np.array([True,False,True,False])) + def test_one_half_space(self): + n = np.array([[-1],[0],[0]]) + x0 = np.array([[0],[0],[0]]) + pts = np.array([[1,-1],[0,0],[0,0]]) + out = half_space_int.half_space_int(n,x0,pts) + assert np.all(out == np.array([True,False])) -if __name__ == '__main__': - test_one_half_space() - test_two_half_spaces() +#------------------------------------------------------------------------------# + def test_two_half_spaces(self): + n = np.array([[-1,0],[0,-1],[0,0]]) + x0 = np.array([[0,0],[0,1],[0,0]]) + pts = np.array([[1,-1,1,0],[2,0,2,0],[0,0,0,0]]) + out = half_space_int.half_space_int(n,x0,pts) + assert np.all(out == np.array([True,False,True,False])) + +#------------------------------------------------------------------------------# From 83c5dc58d4fa3d2c178f847cf2e9c9ad0d01fb4c Mon Sep 17 00:00:00 2001 From: Runar Date: Mon, 6 Feb 2017 09:07:53 +0100 Subject: [PATCH 16/42] Changed name: half_space_int.py -> half_space.py Wanted a more general name to add more functions --- half_space_int.py => half_space.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename half_space_int.py => half_space.py (100%) diff --git a/half_space_int.py b/half_space.py similarity index 100% rename from half_space_int.py rename to half_space.py From 073c10fb196f84ef254c590f4c48591d89519369 Mon Sep 17 00:00:00 2001 From: Runar Date: Mon, 6 Feb 2017 09:07:53 +0100 Subject: [PATCH 17/42] Changed name: half_space_int.py -> half_space.py Wanted a more general name to add more functions --- half_space_int.py => half_space.py | 0 tests/{test_half_space_int.py => test_half_space.py} | 6 +++--- 2 files changed, 3 insertions(+), 3 deletions(-) rename half_space_int.py => half_space.py (100%) rename tests/{test_half_space_int.py => test_half_space.py} (85%) diff --git a/half_space_int.py b/half_space.py similarity index 100% rename from half_space_int.py rename to half_space.py diff --git a/tests/test_half_space_int.py b/tests/test_half_space.py similarity index 85% rename from tests/test_half_space_int.py rename to tests/test_half_space.py index c4a7342663..1acf358024 100644 --- a/tests/test_half_space_int.py +++ b/tests/test_half_space.py @@ -1,7 +1,7 @@ import unittest import numpy as np -from utils import half_space_int +from utils import half_space #------------------------------------------------------------------------------# @@ -14,7 +14,7 @@ def test_one_half_space(self): n = np.array([[-1],[0],[0]]) x0 = np.array([[0],[0],[0]]) pts = np.array([[1,-1],[0,0],[0,0]]) - out = half_space_int.half_space_int(n,x0,pts) + out = half_space.half_space_int(n,x0,pts) assert np.all(out == np.array([True,False])) #------------------------------------------------------------------------------# @@ -23,7 +23,7 @@ def test_two_half_spaces(self): n = np.array([[-1,0],[0,-1],[0,0]]) x0 = np.array([[0,0],[0,1],[0,0]]) pts = np.array([[1,-1,1,0],[2,0,2,0],[0,0,0,0]]) - out = half_space_int.half_space_int(n,x0,pts) + out = half_space.half_space_int(n,x0,pts) assert np.all(out == np.array([True,False,True,False])) #------------------------------------------------------------------------------# From 4b21304bf3b829dac2f3b01e5177f51b0f88eda5 Mon Sep 17 00:00:00 2001 From: Alessio Fumagalli Date: Mon, 6 Feb 2017 15:18:44 +0100 Subject: [PATCH 18/42] add a function to compute a point in the intersection of the halfspaces --- half_space.py | 59 +++++++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 57 insertions(+), 2 deletions(-) diff --git a/half_space.py b/half_space.py index f5ad0f7e59..6d96295d3f 100644 --- a/half_space.py +++ b/half_space.py @@ -1,4 +1,5 @@ import numpy as np +import scipy.optimize as opt def half_space_int(n,x0, pts): @@ -8,7 +9,7 @@ def half_space_int(n,x0, pts): Parameters ---------- n : ndarray - This is the normal vectors of the half planes. The normal + This is the normal vectors of the half planes. The normal vectors is assumed to point out of the half spaces. x0 : ndarray Point on the boundary of the half-spaces. Half space i is given @@ -44,5 +45,59 @@ def half_space_int(n,x0, pts): n = np.repeat( n[:,:,np.newaxis], n_pts,axis=2) for i in range(x0.shape[1]): in_hull += np.sum((pts - x0[:,i,:])*n[:,i,:],axis=0)<=0 - + return in_hull==x0.shape[1] + +#------------------------------------------------------------------------------# + +def half_space_pt(n, x0, pts, recompute=True): + """ + Find an interior point for the halfspaces. + + Parameters + ---------- + n : ndarray + This is the normal vectors of the half planes. The normal + vectors is assumed to coherently for all the half spaces + (inward or outward). + x0 : ndarray + Point on the boundary of the half-spaces. Half space i is given + by all points satisfying (x - x0[:,i])*n[:,i]<=0 + pts : ndarray + Points which defines a bounds for the algorithm. + recompute: bool + If the algorithm fails try again with flipped normals. + + Returns + ------- + out: array + Interior point of the halfspaces. + + We use linear programming to find one interior point for the half spaces. + Assume, n halfspaces defined by: aj*x1+bj*x2+cj*x3+dj<=0, j=1..n. + Perform the following linear program: + max(x5) aj*x1+bj*x2+cj*x3+dj*x4+x5<=0, j=1..n + + Then, if [x1,x2,x3,x4,x5] is an optimal solution with x4>0 and x5>0 we get: + aj*(x1/x4)+bj*(x2/x4)+cj*(x3/x4)+dj<=(-x5/x4) j=1..n and (-x5/x4)<0, + and conclude that the point [x1/x4,x2/x4,x3/x4] is in the interior of all + the halfspaces. Since x5 is optimal, this point is "way in" the interior + (good for precision errors). + http://www.qhull.org/html/qhalf.htm#notes + + """ + dim = (1,n.shape[1]) + c = np.array([0,0,0,0,-1]) + A_ub = np.concatenate( (n, [np.sum(-n*x0, axis=0)], np.ones(dim)) ).T + bounds = ( (np.amin(pts[0,:]), np.amax(pts[0,:]) ), + (np.amin(pts[1,:]), np.amax(pts[1,:]) ), + (np.amin(pts[2,:]), np.amax(pts[2,:]) ), + (None, None), (None, None) ) + res = opt.linprog(c, A_ub, np.zeros(dim).T, bounds=bounds) + if recompute and (res.status != 0 or res.x[3] <=0 or res.x[4] <= 0): + return half_space_pt(-n, x0, pts, False) + else: + assert res.status == 0 and res.x[3] > 0 and res.x[4] > 0 + return np.array(res.x[0:3]) / res.x[3] + +#------------------------------------------------------------------------------# From 9ed0b19c9bb86cc3f949878f6cdf261f789ce633 Mon Sep 17 00:00:00 2001 From: Alessio Fumagalli Date: Mon, 6 Feb 2017 15:19:10 +0100 Subject: [PATCH 19/42] add test to validate the point in the halfspace intersection --- tests/test_half_space.py | 38 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 38 insertions(+) diff --git a/tests/test_half_space.py b/tests/test_half_space.py index 1acf358024..765b144080 100644 --- a/tests/test_half_space.py +++ b/tests/test_half_space.py @@ -26,4 +26,42 @@ def test_two_half_spaces(self): out = half_space.half_space_int(n,x0,pts) assert np.all(out == np.array([True,False,True,False])) +#------------------------------------------------------------------------------# + + def test_half_space_pt_convex_2d(self): + n = np.array([[0,0,0,0,1,-1],[1,1,-1,-1,0,0],[0,0,0,0,0,0]]) + x0 = np.array([[0,0,0,0,0,2/3],[0,0,1/3,1/3,0,0],[0,0,0,0,0,0]]) + pts = np.array([[0,2/3],[0,1/3],[0,0]]) + pt = half_space_int.half_space_pt(n, x0, pts) + assert np.allclose(pt,[1/6,1/6,0]) + +#------------------------------------------------------------------------------# + + def test_half_space_pt_star_shaped_2d(self): + n = np.array([[0,0,1,0,-1,0,1],[1,1,0,1,0,-1,0],[0,0,0,0,0,0,0]]) + x0 = np.array([[0,0,2/3,0,1,0,0],[1/3,1/3,0,0,0,2/3,0],[0,0,0,0,0,0,0]]) + pts = np.array([[0,1],[0,2/3],[0,0]]) + pt = half_space_int.half_space_pt(n, x0, pts) + assert np.allclose(pt,[5/6,1/2,0]) + +#------------------------------------------------------------------------------# + + def test_half_space_pt_convex_3d(self): + n = np.array([[1,-1,0,0,0,0],[0,0,1,-1,0,0],[0,0,0,0,1,-1]]) + x0 = np.array([[0,1,0,0,0,0],[0,0,0,1,0,0],[0,0,0,0,0,1]]) + pts = np.array([[1,0,0],[0,1,0],[0,0,1]]) + pt = half_space_int.half_space_pt(n, x0, pts) + assert np.allclose(pt,[1/2,1/2,1/2]) + +#------------------------------------------------------------------------------# + + def test_half_space_pt_star_shaped_3d(self): + n = np.array([[0,0,1,0,-1,0,1,0,0],[1,1,0,1,0,-1,0,0,0], + [0,0,0,0,0,0,0,1,-1]]) + x0 = np.array([[0,0,2/3,0,1,0,0,0,0],[1/3,1/3,0,0,0,2/3,0,0,0], + [0,0,0,0,0,0,0,0,1]]) + pts = np.array([[0,1],[0,2/3],[0,1]]) + pt = half_space_int.half_space_pt(n, x0, pts) + assert np.allclose(pt,[5/6,1/2,1/6]) + #------------------------------------------------------------------------------# From 27bea32b9cbd2e85b4c526e417fd2474c6622be8 Mon Sep 17 00:00:00 2001 From: Alessio Fumagalli Date: Mon, 6 Feb 2017 15:19:42 +0100 Subject: [PATCH 20/42] add numpy.unique version 1.13.0 (not yet included in anaconda) --- unique.py | 124 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 124 insertions(+) create mode 100644 unique.py diff --git a/unique.py b/unique.py new file mode 100644 index 0000000000..0708ce602d --- /dev/null +++ b/unique.py @@ -0,0 +1,124 @@ +import numpy as np + +def unique_np1130(ar, return_index=False, return_inverse=False, + return_counts=False, axis=None): + """ + Find the unique elements of an array. + Returns the sorted unique elements of an array. There are three optional + outputs in addition to the unique elements: the indices of the input array + that give the unique values, the indices of the unique array that + reconstruct the input array, and the number of times each unique value + comes up in the input array. + Parameters + ---------- + ar : array_like + Input array. Unless `axis` is specified, this will be flattened if it + is not already 1-D. + return_index : bool, optional + If True, also return the indices of `ar` (along the specified axis, + if provided, or in the flattened array) that result in the unique array. + return_inverse : bool, optional + If True, also return the indices of the unique array (for the specified + axis, if provided) that can be used to reconstruct `ar`. + return_counts : bool, optional + If True, also return the number of times each unique item appears + in `ar`. + .. versionadded:: 1.9.0 + axis : int or None, optional + The axis to operate on. If None, `ar` will be flattened beforehand. + Otherwise, duplicate items will be removed along the provided axis, + with all the other axes belonging to the each of the unique elements. + Object arrays or structured arrays that contain objects are not + supported if the `axis` kwarg is used. + .. versionadded:: 1.13.0 + Returns + ------- + unique : ndarray + The sorted unique values. + unique_indices : ndarray, optional + The indices of the first occurrences of the unique values in the + original array. Only provided if `return_index` is True. + unique_inverse : ndarray, optional + The indices to reconstruct the original array from the + unique array. Only provided if `return_inverse` is True. + unique_counts : ndarray, optional + The number of times each of the unique values comes up in the + original array. Only provided if `return_counts` is True. + .. versionadded:: 1.9.0 + See Also + -------- + numpy.lib.arraysetops : Module with a number of other functions for + performing set operations on arrays. + Examples + -------- + >>> np.unique([1, 1, 2, 2, 3, 3]) + array([1, 2, 3]) + >>> a = np.array([[1, 1], [2, 3]]) + >>> np.unique(a) + array([1, 2, 3]) + Return the unique rows of a 2D array + >>> a = np.array([[1, 0, 0], [1, 0, 0], [2, 3, 4]]) + >>> np.unique(a, axis=0) + array([[1, 0, 0], [2, 3, 4]]) + Return the indices of the original array that give the unique values: + >>> a = np.array(['a', 'b', 'b', 'c', 'a']) + >>> u, indices = np.unique(a, return_index=True) + >>> u + array(['a', 'b', 'c'], + dtype='|S1') + >>> indices + array([0, 1, 3]) + >>> a[indices] + array(['a', 'b', 'c'], + dtype='|S1') + Reconstruct the input array from the unique values: + >>> a = np.array([1, 2, 6, 4, 2, 3, 2]) + >>> u, indices = np.unique(a, return_inverse=True) + >>> u + array([1, 2, 3, 4, 6]) + >>> indices + array([0, 1, 4, 3, 1, 2, 1]) + >>> u[indices] + array([1, 2, 6, 4, 2, 3, 2]) + """ + ar = np.asanyarray(ar) + if axis is None: + return np.unique(ar, return_index, return_inverse, return_counts) + if not (-ar.ndim <= axis < ar.ndim): + raise ValueError('Invalid axis kwarg specified for unique') + + ar = np.swapaxes(ar, axis, 0) + orig_shape, orig_dtype = ar.shape, ar.dtype + # Must reshape to a contiguous 2D array for this to work... + ar = ar.reshape(orig_shape[0], -1) + ar = np.ascontiguousarray(ar) + + if ar.dtype.char in (np.typecodes['AllInteger'] + + np.typecodes['Datetime'] + 'S'): + # Optimization: Creating a view of your data with a np.void data type of + # size the number of bytes in a full row. Handles any type where items + # have a unique binary representation, i.e. 0 is only 0, not +0 and -0. + dtype = np.dtype((np.void, ar.dtype.itemsize * ar.shape[1])) + else: + dtype = [('f{i}'.format(i=i), ar.dtype) for i in range(ar.shape[1])] + + try: + consolidated = ar.view(dtype) + except TypeError: + # There's no good way to do this for object arrays, etc... + msg = 'The axis argument to unique is not supported for dtype {dt}' + raise TypeError(msg.format(dt=ar.dtype)) + + def reshape_uniq(uniq): + uniq = uniq.view(orig_dtype) + uniq = uniq.reshape(-1, *orig_shape[1:]) + uniq = np.swapaxes(uniq, 0, axis) + return uniq + + output = np.unique(consolidated, return_index, + return_inverse, return_counts) + if not (return_index or return_inverse or return_counts): + return reshape_uniq(output) + else: + uniq = reshape_uniq(output[0]) + return (uniq,) + output[1:] From 7722a7a37c6bf5be4438eb4c8ce276b935c58d91 Mon Sep 17 00:00:00 2001 From: Alessio Fumagalli Date: Tue, 7 Feb 2017 17:41:46 +0100 Subject: [PATCH 21/42] update name in the test --- tests/test_half_space.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/tests/test_half_space.py b/tests/test_half_space.py index 765b144080..0e96391a9a 100644 --- a/tests/test_half_space.py +++ b/tests/test_half_space.py @@ -32,7 +32,7 @@ def test_half_space_pt_convex_2d(self): n = np.array([[0,0,0,0,1,-1],[1,1,-1,-1,0,0],[0,0,0,0,0,0]]) x0 = np.array([[0,0,0,0,0,2/3],[0,0,1/3,1/3,0,0],[0,0,0,0,0,0]]) pts = np.array([[0,2/3],[0,1/3],[0,0]]) - pt = half_space_int.half_space_pt(n, x0, pts) + pt = half_space.half_space_pt(n, x0, pts) assert np.allclose(pt,[1/6,1/6,0]) #------------------------------------------------------------------------------# @@ -41,7 +41,7 @@ def test_half_space_pt_star_shaped_2d(self): n = np.array([[0,0,1,0,-1,0,1],[1,1,0,1,0,-1,0],[0,0,0,0,0,0,0]]) x0 = np.array([[0,0,2/3,0,1,0,0],[1/3,1/3,0,0,0,2/3,0],[0,0,0,0,0,0,0]]) pts = np.array([[0,1],[0,2/3],[0,0]]) - pt = half_space_int.half_space_pt(n, x0, pts) + pt = half_space.half_space_pt(n, x0, pts) assert np.allclose(pt,[5/6,1/2,0]) #------------------------------------------------------------------------------# @@ -50,7 +50,7 @@ def test_half_space_pt_convex_3d(self): n = np.array([[1,-1,0,0,0,0],[0,0,1,-1,0,0],[0,0,0,0,1,-1]]) x0 = np.array([[0,1,0,0,0,0],[0,0,0,1,0,0],[0,0,0,0,0,1]]) pts = np.array([[1,0,0],[0,1,0],[0,0,1]]) - pt = half_space_int.half_space_pt(n, x0, pts) + pt = half_space.half_space_pt(n, x0, pts) assert np.allclose(pt,[1/2,1/2,1/2]) #------------------------------------------------------------------------------# @@ -61,7 +61,7 @@ def test_half_space_pt_star_shaped_3d(self): x0 = np.array([[0,0,2/3,0,1,0,0,0,0],[1/3,1/3,0,0,0,2/3,0,0,0], [0,0,0,0,0,0,0,0,1]]) pts = np.array([[0,1],[0,2/3],[0,1]]) - pt = half_space_int.half_space_pt(n, x0, pts) + pt = half_space.half_space_pt(n, x0, pts) assert np.allclose(pt,[5/6,1/2,1/6]) #------------------------------------------------------------------------------# From e4c8ef57becae5ee6926d676252eff1a057a37f3 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Wed, 8 Feb 2017 05:03:51 +0100 Subject: [PATCH 22/42] Unit test for setmembership uses standard unittest framework. --- tests/test_setmembership.py | 65 ++++++++++++++++++++----------------- 1 file changed, 36 insertions(+), 29 deletions(-) diff --git a/tests/test_setmembership.py b/tests/test_setmembership.py index 94a9d11709..733a8d4d45 100644 --- a/tests/test_setmembership.py +++ b/tests/test_setmembership.py @@ -5,46 +5,53 @@ @author: keile """ import numpy as np +import unittest from utils import setmembership -def test_unique_rows_1(): - a = np.array([[1, 2], [2, 1], [2, 4], [2, 1], [2, 4]]) - ua_expected = np.array([[1, 2], [2, 1], [2, 4]]) - ia_expected = np.array([0, 1, 2]) - ic_expected = np.array([0, 1, 2, 1, 2]) - ua, ia, ic = setmembership.unique_rows(a) - assert np.sum(np.abs(ua) - np.abs(ua_expected)) == 0 - assert np.all(ia - ia_expected == 0) - assert np.all(ic - ic_expected == 0) +class TestUniqueRows(unittest.TestCase): + def test_unique_rows_1(self): + a = np.array([[1, 2], [2, 1], [2, 4], [2, 1], [2, 4]]) + ua_expected = np.array([[1, 2], [2, 1], [2, 4]]) + ia_expected = np.array([0, 1, 2]) + ic_expected = np.array([0, 1, 2, 1, 2]) + ua, ia, ic = setmembership.unique_rows(a) + assert np.sum(np.abs(ua) - np.abs(ua_expected)) == 0 + assert np.all(ia - ia_expected == 0) + assert np.all(ic - ic_expected == 0) -def test_ismember_rows_with_sort(): - a = np.array([[1, 3, 3, 1, 7], [3, 3, 2, 3, 0]]) - b = np.array([[3, 1, 3, 5, 3], [3, 3, 2, 1, 2]]) - ma, ia = setmembership.ismember_rows(a, b) + if __name__ == '__main__': + unittest.main() - ma_known = np.array([1, 1, 1, 1, 0], dtype=bool) - ia_known = np.array([1, 0, 2, 1]) - assert np.allclose(ma, ma_known) - assert np.allclose(ia, ia_known) +class TestIsmember(unittest.TestCase): + def test_ismember_rows_with_sort(self): + a = np.array([[1, 3, 3, 1, 7], [3, 3, 2, 3, 0]]) + b = np.array([[3, 1, 3, 5, 3], [3, 3, 2, 1, 2]]) + ma, ia = setmembership.ismember_rows(a, b) -def test_ismember_rows_no_sort(): - a = np.array([[1, 3, 3, 1, 7], [3, 3, 2, 3, 0]]) - b = np.array([[3, 1, 2, 5, 3], [3, 3, 3, 1, 1]]) - ma, ia = setmembership.ismember_rows(a, b, sort=False) + ma_known = np.array([1, 1, 1, 1, 0], dtype=bool) + ia_known = np.array([1, 0, 2, 1]) - ma_known = np.array([1, 1, 0, 1, 0], dtype=bool) - ia_known = np.array([1, 0, 1]) + assert np.allclose(ma, ma_known) + assert np.allclose(ia, ia_known) - assert np.allclose(ma, ma_known) - assert np.allclose(ia, ia_known) + + def test_ismember_rows_no_sort(self): + a = np.array([[1, 3, 3, 1, 7], [3, 3, 2, 3, 0]]) + b = np.array([[3, 1, 2, 5, 3], [3, 3, 3, 1, 1]]) + ma, ia = setmembership.ismember_rows(a, b, sort=False) + + ma_known = np.array([1, 1, 0, 1, 0], dtype=bool) + ia_known = np.array([1, 0, 1]) + + assert np.allclose(ma, ma_known) + assert np.allclose(ia, ia_known) + + if __name__ == '__main__': + unittest.main() -if __name__ == '__main__': - test_unique_rows_1() - test_ismember_rows_with_sort() - test_ismember_rows_no_sort() \ No newline at end of file From 0ff4c893a0a21e95d0fcda24f2c8aae34617885b Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Wed, 8 Feb 2017 05:20:15 +0100 Subject: [PATCH 23/42] Function to uniqify matrix columns using tolerance for equality. The function resembles the matlab uniquetol functionality. --- setmembership.py | 79 +++++++++++++++++++++++++++++++++++++ tests/test_setmembership.py | 28 +++++++++++++ 2 files changed, 107 insertions(+) diff --git a/setmembership.py b/setmembership.py index ac15068a32..f0d02283e3 100644 --- a/setmembership.py +++ b/setmembership.py @@ -12,10 +12,14 @@ def unique_rows(data): """ Function similar to Matlab's unique(...,'rows') + See also function unique_columns in this module; this is likely slower, but + is understandable, documented, and has a tolerance option. + Copied from http://stackoverflow.com/questions/16970982/find-unique-rows-in-numpy-array/ (summary pretty far down on the page) Note: I have no idea what happens here + """ b = np.ascontiguousarray(data).view(np.dtype((np.void, data.dtype.itemsize * data.shape[1]))) @@ -102,3 +106,78 @@ def ismember_rows(a, b, sort=True): ind_of_a_in_b = find_occ(occ_a, occ_b) return ismem_a, ind_of_a_in_b + +#--------------------------------------------------------- + +def unique_columns_tol(mat, tol=1e-8, exponent=2): + """ + Remove duplicates from a point set, for a given distance traveling. + + Resembles Matlab's uniquetol function, as applied to columns. To rather + work at rows, use a transpose. + + + Parameters: + mat (np.ndarray, nd x n_pts): Columns to be uniquified + tol (double, optional): Tolerance for when columns are considered equal. + Should be seen in connection with distance between the points in + the points (due to rounding errors). Defaults to 1e-8. + exponent (double, optional): Exponnet in norm used in distance + calculation. Defaults to 2. + + Returns: + np.ndarray: Unique columns. + new_2_old: Index of which points that are preserved + old_2_new: Index of the representation of old points in the reduced + list. + + Example (won't work as doctest): + >>> p_un, n2o, o2n = unique_columns(np.array([[1, 0, 1], [1, 0, 1]])) + >>> p_un + array([[1, 0], [1, 0]]) + >>> n2o + array([0, 1]) + >>> o2n + array([0, 1, 0]) + + """ + + def dist(p, pset): + " Helper function to compute distance " + if p.ndim == 1: + pt = p.reshape((-1, 1)) + else: + pt = p + + return np.power(np.sum(np.power(np.abs(pt - pset), exponent), + axis=0), 1/exponent) + + (nd, l) = mat.shape + + # By default, no columns are kept + keep = np.zeros(l, dtype=np.bool) + + # We will however keep the first point + keep[0] = True + keep_counter = 1 + + # Map from old points to the unique subspace. Defaults to map to itself. + old_2_new = np.arange(l) + + # Loop over all points, check if it is already represented in the kept list + for i in range(1, l): + proximate = np.argwhere(dist(mat[:, i], mat[:, keep]) < tol * np.sqrt(nd)) + + if proximate.size > 0: + # We will not keep this point + old_2_new[i] = proximate[0] + else: + # We have found a new point + keep[i] = True + old_2_new[i] = keep_counter + keep_counter += 1 + # Finally find which elements we kept + new_2_old = np.argwhere(keep).ravel() + + return mat[:, keep], new_2_old, old_2_new + diff --git a/tests/test_setmembership.py b/tests/test_setmembership.py index 733a8d4d45..a35fdf7906 100644 --- a/tests/test_setmembership.py +++ b/tests/test_setmembership.py @@ -55,3 +55,31 @@ def test_ismember_rows_no_sort(self): unittest.main() +class TestUniqueColumns(unittest.TestCase): + + def test_no_common_points(self): + p = np.array([[0, 1, 2], [0, 0, 0]]) + p_unique, new_2_old, old_2_new = setmembership.unique_columns_tol(p) + + assert np.allclose(p, p_unique) + assert np.alltrue(old_2_new == np.arange(3)) + assert np.alltrue(old_2_new == new_2_old) + + def test_remove_one_point(self): + p = np.ones((2, 2)) + p_unique, new_2_old, old_2_new = setmembership.unique_columns_tol(p) + + assert np.allclose(p, np.ones((2, 1))) + assert np.alltrue(old_2_new == np.zeros(2)) + assert np.alltrue(new_2_old == np.zeros(1)) + + def test_remove_one_of_tree(self): + p = np.array([[1, 1, 0], [1, 1, 0]]) + p_unique, new_2_old, old_2_new = setmembership.unique_columns_tol(p) + + assert np.allclose(p_unique, np.array([[1, 0], [1, 0]])) + assert np.alltrue(old_2_new == np.array([0, 0, 1])) + assert np.alltrue(new_2_old == np.array([0, 2])) + + if __name__ == '__main__': + unittest.main() From fe96b2551d363e64a6789d4330e996880568f876 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Wed, 8 Feb 2017 05:34:18 +0100 Subject: [PATCH 24/42] Improved comments in setmembership.ismember_rows. The function should be renamed, return to this later. --- setmembership.py | 48 +++++++++++++++++++++++++++++++++--------------- 1 file changed, 33 insertions(+), 15 deletions(-) diff --git a/setmembership.py b/setmembership.py index f0d02283e3..e5046220f7 100644 --- a/setmembership.py +++ b/setmembership.py @@ -29,7 +29,7 @@ def unique_rows(data): return data[ia], ia, ic -def asvoid(arr): +def _asvoid(arr): """ Taken from @@ -48,7 +48,7 @@ def asvoid(arr): return arr.view(np.dtype((np.void, arr.dtype.itemsize * arr.shape[-1]))) -def find_occ(a, b): +def _find_occ(a, b): """ @@ -73,16 +73,34 @@ def find_occ(a, b): def ismember_rows(a, b, sort=True): """ - Examples - >>> a = np.array([[1, 3, 3, 1, 7], [3, 3, 2, 3, 0]]) - >>> b = np.array([[3, 1, 3, 5, 3], [3, 3, 2, 1, 2]]) - >>> ismember_rows(a, b) - (array([ True, True, True, True, False], dtype=bool), [1, 0, 2, 1]) - - >>> a = np.array([[1, 3, 3, 1, 7], [3, 3, 2, 3, 0]]) - >>> b = np.array([[3, 1, 2, 5, 1], [3, 3, 3, 1, 2]]) - >>> ismember_rows(a, b, sort=False) - (array([ True, True, False, True, False], dtype=bool), [1, 0, 1]) + Find *columns* of a that are also members of *columns* of b. + + The function mimics Matlab's function ismember(..., 'rows'). + + TODO: Rename function, this is confusing! + + Parameters: + a (np.array): Each column in a will search for an equal in b. + b (np.array): Array in which we will look for a twin + sort (boolean, optional): If true, the arrays will be sorted before + seraching, increasing the chances for a match. Defaults to True. + + Returns: + np.array (boolean): For each column in a, true if there is a + corresponding column in b. + np.array (int): Indexes so that b[:, ind] is also found in a. + + Examples: + >>> a = np.array([[1, 3, 3, 1, 7], [3, 3, 2, 3, 0]]) + >>> b = np.array([[3, 1, 3, 5, 3], [3, 3, 2, 1, 2]]) + >>> ismember_rows(a, b) + (array([ True, True, True, True, False], dtype=bool), [1, 0, 2, 1]) + + >>> a = np.array([[1, 3, 3, 1, 7], [3, 3, 2, 3, 0]]) + >>> b = np.array([[3, 1, 2, 5, 1], [3, 3, 3, 1, 2]]) + >>> ismember_rows(a, b, sort=False) + (array([ True, True, False, True, False], dtype=bool), [1, 0, 1]) + """ if sort: sa = np.sort(a, axis=0) @@ -91,8 +109,8 @@ def ismember_rows(a, b, sort=True): sa = a sb = b - voida = asvoid(sa.transpose()) - voidb = asvoid(sb.transpose()) + voida = _asvoid(sa.transpose()) + voidb = _asvoid(sb.transpose()) unq, j, k, count = np.unique(np.vstack((voida, voidb)), return_index=True, return_inverse=True, return_counts=True) @@ -104,7 +122,7 @@ def ismember_rows(a, b, sort=True): occ_a = k[ind_a[ismem_a]] occ_b = k[ind_b] - ind_of_a_in_b = find_occ(occ_a, occ_b) + ind_of_a_in_b = _find_occ(occ_a, occ_b) return ismem_a, ind_of_a_in_b #--------------------------------------------------------- From e45da3f9be7cb80cc7588e442d647d2b966d3cdc Mon Sep 17 00:00:00 2001 From: Runar Date: Fri, 10 Feb 2017 13:55:37 +0100 Subject: [PATCH 25/42] Added graph class with unit tests The graph class create a graph and is able to color the nodes int the graph by the connected regions. --- graph.py | 59 +++++++++++++++++++++++++++++++++++++++++++++ tests/test_graph.py | 57 +++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 116 insertions(+) create mode 100644 graph.py create mode 100644 tests/test_graph.py diff --git a/graph.py b/graph.py new file mode 100644 index 0000000000..d793ef3150 --- /dev/null +++ b/graph.py @@ -0,0 +1,59 @@ +import numpy as np +from scipy import sparse as sps + +class Graph: + """ + A graph class. + + The graph class stores the nodes and edges of the graph in a sparse + array (equivalently to face_nodes in the Grid class). + + Attributes: + node_connections (sps.csc-matrix): Should be given at construction. + node_node connections. Matrix size: num_nodes x num_nodes. + node_connections[i,j] should be true if there is an edge + connecting node i and j. + regions (int) the number of regions. A region is a set of nodes + that can be reached by traversing the graph. Two nodes are + int different regions if they can not be reached by traversing + the graph. + color (int) the color of each region. Initialized as (NaN). By + calling color_nodes() all nodes in a region are given the + same colors and nodes in different regions are given different + colors. + """ + def __init__(self,node_connections): + self.node_connections = node_connections + self.regions = 0 + self.color = np.array([np.NaN]*node_connections.shape[1]) + + def color_nodes(self): + """ + Color the nodes in each region + """ + color = 0 + while self.regions Date: Sun, 19 Feb 2017 14:13:19 +0100 Subject: [PATCH 26/42] add functions to compute normL2 and errorL2 --- errors/error.py | 106 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 106 insertions(+) create mode 100644 errors/error.py diff --git a/errors/error.py b/errors/error.py new file mode 100644 index 0000000000..df8b5d8540 --- /dev/null +++ b/errors/error.py @@ -0,0 +1,106 @@ +import numpy as np + +#------------------------------------------------------------------------------# + +def interpolate(g, fun): + """ + Interpolate a scalar or vector function on the cell centers of the grid. + + Parameters + ---------- + g : grid + Grid, or a subclass, with geometry fields computed. + fun : function + Scalar or vector function. + + Return + ------ + out: np.ndarray (dim of fun, g.num_cells) + Function interpolated in the cell centers. + + Examples + -------- + def fun_p(pt): return np.sin(2*np.pi*pt[0])*np.sin(2*np.pi*pt[1]) + def fun_u(pt): return [\ + -2*np.pi*np.cos(2*np.pi*pt[0])*np.sin(2*np.pi*pt[1]), + -2*np.pi*np.sin(2*np.pi*pt[0])*np.cos(2*np.pi*pt[1])] + p_ex = interpolate(g, fun_p) + u_ex = interpolate(g, fun_u) + + """ + + return np.array([fun(pt) for pt in g.cell_centers.T]).T + +#------------------------------------------------------------------------------# + +def norm_L2(g, val): + """ + Compute the L2 norm of a scalar or vector field. + + Parameters + ---------- + g : grid + Grid, or a subclass, with geometry fields computed. + val : np.ndarray (dim of val, g.num_cells) + Scalar or vector field. + + Return + ------ + out: double + The L2 norm of the input field. + + Examples + -------- + def fun_p(pt): return np.sin(2*np.pi*pt[0])*np.sin(2*np.pi*pt[1]) + p_ex = interpolate(g, fun_p) + norm_ex = norm_L2(g, p_ex) + + """ + + val = np.asarray(val) + norm_sq = lambda v: np.sum(np.multiply(np.square(v), g.cell_volumes)) + if val.ndim == 1: + return np.sqrt(norm_sq(val)) + else: + return np.sqrt(np.sum([norm_sq(v) for v in val])) + +#------------------------------------------------------------------------------# + +def error_L2(g, val, val_ex, relative=True): + """ + Compute the L2 error of a scalar or vector field with respect to a reference + field. It is possible to compute the relative error (default) or the + absolute error. + + Parameters + ---------- + g : grid + Grid, or a subclass, with geometry fields computed. + val : np.ndarray (dim of val, g.num_cells) + Scalar or vector field. + val_ex : np.ndarray (dim of val, g.num_cells) + Reference scalar or vector field, i.e. the exact solution + relative: bool (True default) + Compute the relative error (if True) or the absolute error (if False). + + Return + ------ + out: double + The L2 error of the input fields. + + Examples + -------- + p = ... + def fun_p(pt): return np.sin(2*np.pi*pt[0])*np.sin(2*np.pi*pt[1]) + p_ex = interpolate(g, fun_p) + err_p = err_L2(g, p, p_ex) + + """ + + val, val_ex = np.asarray(val), np.asarray(val_ex) + err = norm_L2(g, np.subtract(val, val_ex)) + den = norm_L2(g, val_ex) if relative else 1 + assert den != 0 + return err/den + +#------------------------------------------------------------------------------# From 3039f5f74f24206c8a9cdf51cd761a87c8c51911 Mon Sep 17 00:00:00 2001 From: Alessio Fumagalli Date: Wed, 22 Feb 2017 13:50:07 +0100 Subject: [PATCH 27/42] fix indentation --- graph.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/graph.py b/graph.py index d793ef3150..8c107f5e8d 100644 --- a/graph.py +++ b/graph.py @@ -7,19 +7,19 @@ class Graph: The graph class stores the nodes and edges of the graph in a sparse array (equivalently to face_nodes in the Grid class). - + Attributes: node_connections (sps.csc-matrix): Should be given at construction. - node_node connections. Matrix size: num_nodes x num_nodes. - node_connections[i,j] should be true if there is an edge + node_node connections. Matrix size: num_nodes x num_nodes. + node_connections[i,j] should be true if there is an edge connecting node i and j. - regions (int) the number of regions. A region is a set of nodes - that can be reached by traversing the graph. Two nodes are + regions (int) the number of regions. A region is a set of nodes + that can be reached by traversing the graph. Two nodes are int different regions if they can not be reached by traversing the graph. - color (int) the color of each region. Initialized as (NaN). By + color (int) the color of each region. Initialized as (NaN). By calling color_nodes() all nodes in a region are given the - same colors and nodes in different regions are given different + same colors and nodes in different regions are given different colors. """ def __init__(self,node_connections): From 09e8be8c7c823a73164905d44cb8e10abb5033c9 Mon Sep 17 00:00:00 2001 From: Alessio Fumagalli Date: Wed, 22 Feb 2017 13:50:14 +0100 Subject: [PATCH 28/42] remove print --- graph.py | 1 - 1 file changed, 1 deletion(-) diff --git a/graph.py b/graph.py index 8c107f5e8d..d2adcc72ed 100644 --- a/graph.py +++ b/graph.py @@ -34,7 +34,6 @@ def color_nodes(self): color = 0 while self.regions Date: Wed, 1 Mar 2017 11:42:45 +0100 Subject: [PATCH 29/42] Bugfix in setmembership --- setmembership.py | 425 +++++++++++++++++++++++++---------------------- 1 file changed, 224 insertions(+), 201 deletions(-) diff --git a/setmembership.py b/setmembership.py index e5046220f7..79c3282166 100644 --- a/setmembership.py +++ b/setmembership.py @@ -1,201 +1,224 @@ -# -*- coding: utf-8 -*- -""" -Created on Thu Feb 25 20:16:42 2016 - -@author: keile -""" - -import numpy as np - - -def unique_rows(data): - """ - Function similar to Matlab's unique(...,'rows') - - See also function unique_columns in this module; this is likely slower, but - is understandable, documented, and has a tolerance option. - - Copied from - http://stackoverflow.com/questions/16970982/find-unique-rows-in-numpy-array/ - (summary pretty far down on the page) - Note: I have no idea what happens here - - """ - b = np.ascontiguousarray(data).view(np.dtype((np.void, - data.dtype.itemsize * data.shape[1]))) - _, ia = np.unique(b, return_index=True) - _, ic = np.unique(b, return_inverse=True) -# return np.unique(b).view(data.dtype).reshape(-1, data.shape[1]), ia, ic - return data[ia], ia, ic - - -def _asvoid(arr): - """ - - Taken from - http://stackoverflow.com/questions/22699756/python-version-of-ismember-with-rows-and-index - - View the array as dtype np.void (bytes) - This views the last axis of ND-arrays as bytes so you can perform - comparisons on the entire row. - http://stackoverflow.com/a/16840350/190597 (Jaime, 2013-05) - Warning: When using asvoid for comparison, note that float zeros may - compare UNEQUALLY - >>> asvoid([-0.]) == asvoid([0.]) - array([False], dtype=bool) - """ - arr = np.ascontiguousarray(arr) - return arr.view(np.dtype((np.void, arr.dtype.itemsize * arr.shape[-1]))) - - -def _find_occ(a, b): - """ - - - Code snippet found at - http://stackoverflow.com/questions/15864082/python-equivalent-of-matlabs-ismember-function?rq=1 - - Parameters - ---------- - a - b - - Returns - ------- - - """ - bind = {} - for i, elt in enumerate(b): - if elt not in bind: - bind[elt] = i - return [bind.get(itm, None) for itm in a] - - -def ismember_rows(a, b, sort=True): - """ - Find *columns* of a that are also members of *columns* of b. - - The function mimics Matlab's function ismember(..., 'rows'). - - TODO: Rename function, this is confusing! - - Parameters: - a (np.array): Each column in a will search for an equal in b. - b (np.array): Array in which we will look for a twin - sort (boolean, optional): If true, the arrays will be sorted before - seraching, increasing the chances for a match. Defaults to True. - - Returns: - np.array (boolean): For each column in a, true if there is a - corresponding column in b. - np.array (int): Indexes so that b[:, ind] is also found in a. - - Examples: - >>> a = np.array([[1, 3, 3, 1, 7], [3, 3, 2, 3, 0]]) - >>> b = np.array([[3, 1, 3, 5, 3], [3, 3, 2, 1, 2]]) - >>> ismember_rows(a, b) - (array([ True, True, True, True, False], dtype=bool), [1, 0, 2, 1]) - - >>> a = np.array([[1, 3, 3, 1, 7], [3, 3, 2, 3, 0]]) - >>> b = np.array([[3, 1, 2, 5, 1], [3, 3, 3, 1, 2]]) - >>> ismember_rows(a, b, sort=False) - (array([ True, True, False, True, False], dtype=bool), [1, 0, 1]) - - """ - if sort: - sa = np.sort(a, axis=0) - sb = np.sort(b, axis=0) - else: - sa = a - sb = b - - voida = _asvoid(sa.transpose()) - voidb = _asvoid(sb.transpose()) - unq, j, k, count = np.unique(np.vstack((voida, voidb)), return_index=True, - return_inverse=True, return_counts=True) - - num_a = a.shape[1] - ind_a = np.arange(num_a) - ind_b = num_a + np.arange(b.shape[1]) - num_occ_a = count[k[ind_a]] - ismem_a = (num_occ_a > 1) - occ_a = k[ind_a[ismem_a]] - occ_b = k[ind_b] - - ind_of_a_in_b = _find_occ(occ_a, occ_b) - return ismem_a, ind_of_a_in_b - -#--------------------------------------------------------- - -def unique_columns_tol(mat, tol=1e-8, exponent=2): - """ - Remove duplicates from a point set, for a given distance traveling. - - Resembles Matlab's uniquetol function, as applied to columns. To rather - work at rows, use a transpose. - - - Parameters: - mat (np.ndarray, nd x n_pts): Columns to be uniquified - tol (double, optional): Tolerance for when columns are considered equal. - Should be seen in connection with distance between the points in - the points (due to rounding errors). Defaults to 1e-8. - exponent (double, optional): Exponnet in norm used in distance - calculation. Defaults to 2. - - Returns: - np.ndarray: Unique columns. - new_2_old: Index of which points that are preserved - old_2_new: Index of the representation of old points in the reduced - list. - - Example (won't work as doctest): - >>> p_un, n2o, o2n = unique_columns(np.array([[1, 0, 1], [1, 0, 1]])) - >>> p_un - array([[1, 0], [1, 0]]) - >>> n2o - array([0, 1]) - >>> o2n - array([0, 1, 0]) - - """ - - def dist(p, pset): - " Helper function to compute distance " - if p.ndim == 1: - pt = p.reshape((-1, 1)) - else: - pt = p - - return np.power(np.sum(np.power(np.abs(pt - pset), exponent), - axis=0), 1/exponent) - - (nd, l) = mat.shape - - # By default, no columns are kept - keep = np.zeros(l, dtype=np.bool) - - # We will however keep the first point - keep[0] = True - keep_counter = 1 - - # Map from old points to the unique subspace. Defaults to map to itself. - old_2_new = np.arange(l) - - # Loop over all points, check if it is already represented in the kept list - for i in range(1, l): - proximate = np.argwhere(dist(mat[:, i], mat[:, keep]) < tol * np.sqrt(nd)) - - if proximate.size > 0: - # We will not keep this point - old_2_new[i] = proximate[0] - else: - # We have found a new point - keep[i] = True - old_2_new[i] = keep_counter - keep_counter += 1 - # Finally find which elements we kept - new_2_old = np.argwhere(keep).ravel() - - return mat[:, keep], new_2_old, old_2_new - +# -*- coding: utf-8 -*- +""" +Created on Thu Feb 25 20:16:42 2016 + +@author: keile +""" + +import numpy as np + + +def unique_rows(data): + """ + Function similar to Matlab's unique(...,'rows') + + See also function unique_columns in this module; this is likely slower, but + is understandable, documented, and has a tolerance option. + + Copied from + http://stackoverflow.com/questions/16970982/find-unique-rows-in-numpy-array/ + (summary pretty far down on the page) + Note: I have no idea what happens here + + """ + b = np.ascontiguousarray(data).view(np.dtype((np.void, + data.dtype.itemsize * data.shape[1]))) + _, ia = np.unique(b, return_index=True) + _, ic = np.unique(b, return_inverse=True) +# return np.unique(b).view(data.dtype).reshape(-1, data.shape[1]), ia, ic + return data[ia], ia, ic + + +def _asvoid(arr): + """ + + Taken from + http://stackoverflow.com/questions/22699756/python-version-of-ismember-with-rows-and-index + + View the array as dtype np.void (bytes) + This views the last axis of ND-arrays as bytes so you can perform + comparisons on the entire row. + http://stackoverflow.com/a/16840350/190597 (Jaime, 2013-05) + Warning: When using asvoid for comparison, note that float zeros may + compare UNEQUALLY + >>> asvoid([-0.]) == asvoid([0.]) + array([False], dtype=bool) + """ + arr = np.ascontiguousarray(arr) + return arr.view(np.dtype((np.void, arr.dtype.itemsize * arr.shape[-1]))) + + +def _find_occ(a, b): + """ + + + Code snippet found at + http://stackoverflow.com/questions/15864082/python-equivalent-of-matlabs-ismember-function?rq=1 + + Parameters + ---------- + a + b + + Returns + ------- + + """ + bind = {} + for i, elt in enumerate(b): + if elt not in bind: + bind[elt] = i + return [bind.get(itm, None) for itm in a] + + +def ismember_rows(a, b, sort=True): + """ + Find *columns* of a that are also members of *columns* of b. + + The function mimics Matlab's function ismember(..., 'rows'). + + TODO: Rename function, this is confusing! + + Parameters: + a (np.array): Each column in a will search for an equal in b. + b (np.array): Array in which we will look for a twin + sort (boolean, optional): If true, the arrays will be sorted before + seraching, increasing the chances for a match. Defaults to True. + + Returns: + np.array (boolean): For each column in a, true if there is a + corresponding column in b. + np.array (int): Indexes so that b[:, ind] is also found in a. + + Examples: + >>> a = np.array([[1, 3, 3, 1, 7], [3, 3, 2, 3, 0]]) + >>> b = np.array([[3, 1, 3, 5, 3], [3, 3, 2, 1, 2]]) + >>> ismember_rows(a, b) + (array([ True, True, True, True, False], dtype=bool), [1, 0, 2, 1]) + + >>> a = np.array([[1, 3, 3, 1, 7], [3, 3, 2, 3, 0]]) + >>> b = np.array([[3, 1, 2, 5, 1], [3, 3, 3, 1, 2]]) + >>> ismember_rows(a, b, sort=False) + (array([ True, True, False, True, False], dtype=bool), [1, 0, 1]) + + """ + if sort: + sa = np.sort(a, axis=0) + sb = np.sort(b, axis=0) + else: + sa = a + sb = b + + if a.ndim == 1: + num_a = 1 + else: + num_a = a.shape[1] + + """ + Old code, found on the internet. Not sure about the reliability, so we kick + it out for the moment. + voida = _asvoid(sa.transpose()) + voidb = _asvoid(sb.transpose()) + unq, j, k, count = np.unique(np.vstack((voida, voidb)), return_index=True, + return_inverse=True, return_counts=True) + + ind_a = np.arange(num_a) + ind_b = num_a + np.arange(b.shape[1]) + num_occ_a = count[k[ind_a]] + ismem_a = (num_occ_a > 1) + occ_a = k[ind_a[ismem_a]] + occ_b = k[ind_b] + + ind_of_a_in_b = _find_occ(occ_a, occ_b) + """ + ismem_a = np.zeros(num_a, dtype=np.bool) + ind_of_a_in_b = [] + for i in range(num_a): + if sa.ndim == 1: + diff = np.abs(sb - sa[i]) + else: + diff = np.sum(np.abs(sb - sa[:, i].reshape((-1, 1))), axis=0) + if np.any(diff == 0): + ismem_a[i] = True + ind_of_a_in_b.append(np.where(diff == 0)[0]) + + + # Remove None-types from output indices +# ind_of_a_in_b = [i for i in ind_of_a_in_b if i is not None] + + return ismem_a, np.squeeze(np.array(ind_of_a_in_b, dtype='int')) + +#--------------------------------------------------------- + +def unique_columns_tol(mat, tol=1e-8, exponent=2): + """ + Remove duplicates from a point set, for a given distance traveling. + + Resembles Matlab's uniquetol function, as applied to columns. To rather + work at rows, use a transpose. + + + Parameters: + mat (np.ndarray, nd x n_pts): Columns to be uniquified + tol (double, optional): Tolerance for when columns are considered equal. + Should be seen in connection with distance between the points in + the points (due to rounding errors). Defaults to 1e-8. + exponent (double, optional): Exponnet in norm used in distance + calculation. Defaults to 2. + + Returns: + np.ndarray: Unique columns. + new_2_old: Index of which points that are preserved + old_2_new: Index of the representation of old points in the reduced + list. + + Example (won't work as doctest): + >>> p_un, n2o, o2n = unique_columns(np.array([[1, 0, 1], [1, 0, 1]])) + >>> p_un + array([[1, 0], [1, 0]]) + >>> n2o + array([0, 1]) + >>> o2n + array([0, 1, 0]) + + """ + + def dist(p, pset): + " Helper function to compute distance " + if p.ndim == 1: + pt = p.reshape((-1, 1)) + else: + pt = p + + return np.power(np.sum(np.power(np.abs(pt - pset), exponent), + axis=0), 1/exponent) + + (nd, l) = mat.shape + + # By default, no columns are kept + keep = np.zeros(l, dtype=np.bool) + + # We will however keep the first point + keep[0] = True + keep_counter = 1 + + # Map from old points to the unique subspace. Defaults to map to itself. + old_2_new = np.arange(l) + + # Loop over all points, check if it is already represented in the kept list + for i in range(1, l): + proximate = np.argwhere(dist(mat[:, i], mat[:, keep]) < tol * np.sqrt(nd)) + + if proximate.size > 0: + # We will not keep this point + old_2_new[i] = proximate[0] + else: + # We have found a new point + keep[i] = True + old_2_new[i] = keep_counter + keep_counter += 1 + # Finally find which elements we kept + new_2_old = np.argwhere(keep).ravel() + + return mat[:, keep], new_2_old, old_2_new + From 69f78c08dde82a351045a02b12942fb13e43ede7 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Wed, 1 Mar 2017 14:10:05 +0100 Subject: [PATCH 30/42] Bugfix in setmembership.ismember_rows --- setmembership.py | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/setmembership.py b/setmembership.py index 79c3282166..900a0b5971 100644 --- a/setmembership.py +++ b/setmembership.py @@ -132,7 +132,7 @@ def ismember_rows(a, b, sort=True): ind_of_a_in_b = _find_occ(occ_a, occ_b) """ ismem_a = np.zeros(num_a, dtype=np.bool) - ind_of_a_in_b = [] + ind_of_a_in_b = np.empty(0) for i in range(num_a): if sa.ndim == 1: diff = np.abs(sb - sa[i]) @@ -140,13 +140,12 @@ def ismember_rows(a, b, sort=True): diff = np.sum(np.abs(sb - sa[:, i].reshape((-1, 1))), axis=0) if np.any(diff == 0): ismem_a[i] = True - ind_of_a_in_b.append(np.where(diff == 0)[0]) + hit = np.where(diff==0)[0] + if hit.size > 1: + hit = hit[0] + ind_of_a_in_b = np.append(ind_of_a_in_b, hit) - - # Remove None-types from output indices -# ind_of_a_in_b = [i for i in ind_of_a_in_b if i is not None] - - return ismem_a, np.squeeze(np.array(ind_of_a_in_b, dtype='int')) + return ismem_a, ind_of_a_in_b.astype('int') #--------------------------------------------------------- From 7fa300d57051bb96ecf142791b473fb11025947d Mon Sep 17 00:00:00 2001 From: Runar Date: Thu, 2 Mar 2017 14:02:32 +0100 Subject: [PATCH 31/42] Removed a print statement from graph.py --- graph.py | 1 - 1 file changed, 1 deletion(-) diff --git a/graph.py b/graph.py index d793ef3150..36d20eb88f 100644 --- a/graph.py +++ b/graph.py @@ -34,7 +34,6 @@ def color_nodes(self): color = 0 while self.regions Date: Mon, 6 Mar 2017 10:29:45 +0100 Subject: [PATCH 32/42] Improved comments in utility function in setmembership --- setmembership.py | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/setmembership.py b/setmembership.py index 900a0b5971..8b779285df 100644 --- a/setmembership.py +++ b/setmembership.py @@ -50,24 +50,25 @@ def _asvoid(arr): def _find_occ(a, b): """ + Find index of occurences of a in b. + The function has only been tested on np.arrays, but it should be fairly + general (only require iterables?) Code snippet found at http://stackoverflow.com/questions/15864082/python-equivalent-of-matlabs-ismember-function?rq=1 - Parameters - ---------- - a - b - - Returns - ------- - """ + # Base search on a dictionary bind = {} + # Invert dictionary to create a map from an item in b to the *first* + # occurence of the item. + # NOTE: If we ever need to give the option of returning last index, it + # should require no more than bypassing the if statement. for i, elt in enumerate(b): if elt not in bind: bind[elt] = i + # Use inverse mapping to obtain return [bind.get(itm, None) for itm in a] From c01926e985964d84c7a44870580da7a9ca3d47e3 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Mon, 6 Mar 2017 10:30:06 +0100 Subject: [PATCH 33/42] Bugfix in ismember_rows in setmebership The previous version sometimes returned True based on items in a. Also improved comments. --- setmembership.py | 93 +++++++++++++++++++++++++++++++----------------- 1 file changed, 60 insertions(+), 33 deletions(-) diff --git a/setmembership.py b/setmembership.py index 8b779285df..c84bf0156c 100644 --- a/setmembership.py +++ b/setmembership.py @@ -72,7 +72,7 @@ def _find_occ(a, b): return [bind.get(itm, None) for itm in a] -def ismember_rows(a, b, sort=True): +def ismember_rows(a, b, sort=True, simple_version=False): """ Find *columns* of a that are also members of *columns* of b. @@ -85,6 +85,9 @@ def ismember_rows(a, b, sort=True): b (np.array): Array in which we will look for a twin sort (boolean, optional): If true, the arrays will be sorted before seraching, increasing the chances for a match. Defaults to True. + simple_verion (boolean, optional): Use an alternative implementation + based on a global for loop. The code is slow for large arrays, but + easy to understand. Defaults to False. Returns: np.array (boolean): For each column in a, true if there is a @@ -115,38 +118,62 @@ def ismember_rows(a, b, sort=True): else: num_a = a.shape[1] - """ - Old code, found on the internet. Not sure about the reliability, so we kick - it out for the moment. - voida = _asvoid(sa.transpose()) - voidb = _asvoid(sb.transpose()) - unq, j, k, count = np.unique(np.vstack((voida, voidb)), return_index=True, - return_inverse=True, return_counts=True) - - ind_a = np.arange(num_a) - ind_b = num_a + np.arange(b.shape[1]) - num_occ_a = count[k[ind_a]] - ismem_a = (num_occ_a > 1) - occ_a = k[ind_a[ismem_a]] - occ_b = k[ind_b] - - ind_of_a_in_b = _find_occ(occ_a, occ_b) - """ - ismem_a = np.zeros(num_a, dtype=np.bool) - ind_of_a_in_b = np.empty(0) - for i in range(num_a): - if sa.ndim == 1: - diff = np.abs(sb - sa[i]) - else: - diff = np.sum(np.abs(sb - sa[:, i].reshape((-1, 1))), axis=0) - if np.any(diff == 0): - ismem_a[i] = True - hit = np.where(diff==0)[0] - if hit.size > 1: - hit = hit[0] - ind_of_a_in_b = np.append(ind_of_a_in_b, hit) - - return ismem_a, ind_of_a_in_b.astype('int') + if simple_version: + # Use straightforward search, based on a for loop. This is slow for + # large arrays, but as the alternative implementation is opaque, and + # there has been some doubts on its reliability, this version is kept + # as a safeguard. + ismem_a = np.zeros(num_a, dtype=np.bool) + ind_of_a_in_b = np.empty(0) + for i in range(num_a): + if sa.ndim == 1: + diff = np.abs(sb - sa[i]) + else: + diff = np.sum(np.abs(sb - sa[:, i].reshape((-1, 1))), axis=0) + if np.any(diff == 0): + ismem_a[i] = True + hit = np.where(diff==0)[0] + if hit.size > 1: + hit = hit[0] + ind_of_a_in_b = np.append(ind_of_a_in_b, hit) + + return ismem_a, ind_of_a_in_b.astype('int') + + else: + + # Represent the arrays as voids to facilitate quick comparison + voida = _asvoid(sa.transpose()) + voidb = _asvoid(sb.transpose()) + + # Use unique to count the number of occurences in a + unq, j, k, count = np.unique(np.vstack((voida, voidb)), return_index=True, + return_inverse=True, return_counts=True) + # Also count the number of occurences in voida + _, j_a, k_a, count_a = np.unique(voida, return_index=True, + return_inverse=True, return_counts=True) + + # Index of a and b elements in the combined array + ind_a = np.arange(num_a) + ind_b = num_a + np.arange(b.shape[1]) + + # Count number of occurences in combine array, and in a only + num_occ_a_and_b = count[k[ind_a]] + num_occ_a = count_a[k_a[ind_a]] + + # Subtraction gives number of a in b + num_occ_a_in_b = num_occ_a_and_b - num_occ_a + ismem_a = (num_occ_a_in_b > 0) + + # To get the indices of common elements in a and b, compare the + # elements in k (pointers to elements in the unique combined arary) + occ_a = k[ind_a[ismem_a]] + occ_b = k[ind_b] + + ind_of_a_in_b = _find_occ(occ_a, occ_b) + # Remove None types when no hit was found + ind_of_a_in_b = [i for i in ind_of_a_in_b if i is not None] + + return ismem_a, np.array(ind_of_a_in_b, dtype='int') #--------------------------------------------------------- From 1356e09ba939c1c113b24ec64eff726ea4db9059 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Mon, 6 Mar 2017 10:32:02 +0100 Subject: [PATCH 34/42] Improved unit tests for setmembership.ismember_rows Also some formating. --- tests/test_setmembership.py | 71 ++++++++++++++++++++++++++++++++++--- 1 file changed, 66 insertions(+), 5 deletions(-) diff --git a/tests/test_setmembership.py b/tests/test_setmembership.py index a35fdf7906..e0eb646822 100644 --- a/tests/test_setmembership.py +++ b/tests/test_setmembership.py @@ -29,8 +29,10 @@ def test_unique_rows_1(self): class TestIsmember(unittest.TestCase): def test_ismember_rows_with_sort(self): - a = np.array([[1, 3, 3, 1, 7], [3, 3, 2, 3, 0]]) - b = np.array([[3, 1, 3, 5, 3], [3, 3, 2, 1, 2]]) + a = np.array([[1, 3, 3, 1, 7], + [3, 3, 2, 3, 0]]) + b = np.array([[3, 1, 3, 5, 3], + [3, 3, 2, 1, 2]]) ma, ia = setmembership.ismember_rows(a, b) ma_known = np.array([1, 1, 1, 1, 0], dtype=bool) @@ -41,8 +43,10 @@ def test_ismember_rows_with_sort(self): def test_ismember_rows_no_sort(self): - a = np.array([[1, 3, 3, 1, 7], [3, 3, 2, 3, 0]]) - b = np.array([[3, 1, 2, 5, 3], [3, 3, 3, 1, 1]]) + a = np.array([[1, 3, 3, 1, 7], + [3, 3, 2, 3, 0]]) + b = np.array([[3, 1, 2, 5, 3], + [3, 3, 3, 1, 1]]) ma, ia = setmembership.ismember_rows(a, b, sort=False) ma_known = np.array([1, 1, 0, 1, 0], dtype=bool) @@ -50,7 +54,64 @@ def test_ismember_rows_no_sort(self): assert np.allclose(ma, ma_known) assert np.allclose(ia, ia_known) - + + def test_ismember_rows_unqual_sizes_1(self): + # a larger than b + a = np.array([[1, 3, 3, 1, 7], + [3, 3, 2, 3, 0]]) + b = np.array([[3, 1, 2, 5], + [3, 3, 3, 1]]) + ma, ia = setmembership.ismember_rows(a, b) + + ma_known = np.array([1, 1, 1, 1, 0], dtype=bool) + ia_known = np.array([1, 0, 2, 1]) + + assert np.allclose(ma, ma_known) + assert np.allclose(ia, ia_known) + + def test_ismember_rows_unqual_sizes_1(self): + # b larger than b + a = np.array([[1, 3, 3, 1, 7], + [3, 3, 2, 3, 0]]) + b = np.array([[3, 1, 2, 5, 3, 4, 7], + [3, 3, 3, 1, 9, 9, 9]]) + ma, ia = setmembership.ismember_rows(a, b) + + ma_known = np.array([1, 1, 1, 1, 0], dtype=bool) + ia_known = np.array([1, 0, 2, 1]) + + assert np.allclose(ma, ma_known) + assert np.allclose(ia, ia_known) + + def test_ismember_rows_double_occurence_a_no_b(self): + # There are duplicate occurences in a that are not found in b + a = np.array([[1, 3, 3, 1, 7], + [3, 3, 2, 3, 0]]) + b = np.array([[3, 2, 5], + [3, 3, 1]]) + ma, ia = setmembership.ismember_rows(a, b) + + ma_known = np.array([0, 1, 1, 0, 0], dtype=bool) + ia_known = np.array([0, 1]) + + assert np.allclose(ma, ma_known) + assert np.allclose(ia, ia_known) + + def test_ismember_rows_double_occurence_a_and_b(self): + # There are duplicate occurences in a, and the same item is found in b + a = np.array([[1, 3, 3, 1, 7], + [3, 3, 2, 3, 0]]) + b = np.array([[3, 1, 2, 5, 3], + [3, 3, 3, 1, 1]]) + ma, ia = setmembership.ismember_rows(a, b) + + ma_known = np.array([1, 1, 1, 1, 0], dtype=bool) + ia_known = np.array([1, 0, 2, 1]) + + assert np.allclose(ma, ma_known) + assert np.allclose(ia, ia_known) + + if __name__ == '__main__': unittest.main() From 18db61814a961dd087e84e154eeb9952b39610f1 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Mon, 6 Mar 2017 11:41:05 +0100 Subject: [PATCH 35/42] Fix of setmembership.ismember_rows for 1d arrays --- setmembership.py | 44 +++++++++++++++++++++++-------------- tests/test_setmembership.py | 23 +++++++++++++++++++ 2 files changed, 50 insertions(+), 17 deletions(-) diff --git a/setmembership.py b/setmembership.py index c84bf0156c..31ffcf3276 100644 --- a/setmembership.py +++ b/setmembership.py @@ -60,6 +60,7 @@ def _find_occ(a, b): """ # Base search on a dictionary + bind = {} # Invert dictionary to create a map from an item in b to the *first* # occurence of the item. @@ -106,17 +107,16 @@ def ismember_rows(a, b, sort=True, simple_version=False): (array([ True, True, False, True, False], dtype=bool), [1, 0, 1]) """ - if sort: + + # Sort if required, but not if the input is 1d + if sort and a.ndim > 1: sa = np.sort(a, axis=0) sb = np.sort(b, axis=0) else: sa = a sb = b - if a.ndim == 1: - num_a = 1 - else: - num_a = a.shape[1] + num_a = a.shape[-1] if simple_version: # Use straightforward search, based on a for loop. This is slow for @@ -140,21 +140,31 @@ def ismember_rows(a, b, sort=True, simple_version=False): return ismem_a, ind_of_a_in_b.astype('int') else: - - # Represent the arrays as voids to facilitate quick comparison - voida = _asvoid(sa.transpose()) - voidb = _asvoid(sb.transpose()) - - # Use unique to count the number of occurences in a - unq, j, k, count = np.unique(np.vstack((voida, voidb)), return_index=True, - return_inverse=True, return_counts=True) - # Also count the number of occurences in voida - _, j_a, k_a, count_a = np.unique(voida, return_index=True, - return_inverse=True, return_counts=True) + if a.ndim == 1: + # Special treatment of 1d, vstack of voids (below) ran into trouble + # here. + unq, k, count = np.unique(np.hstack((a, b)), return_inverse=True, + return_counts=True) + _, k_a, count_a = np.unique(a, return_inverse=True, + return_counts=True) + else: + # Represent the arrays as voids to facilitate quick comparison + voida = _asvoid(sa.transpose()) + voidb = _asvoid(sb.transpose()) + + # Use unique to count the number of occurences in a + unq, j, k, count = np.unique(np.vstack((voida, voidb)), + return_index=True, + return_inverse=True, + return_counts=True) + # Also count the number of occurences in voida + _, j_a, k_a, count_a = np.unique(voida, return_index=True, + return_inverse=True, + return_counts=True) # Index of a and b elements in the combined array ind_a = np.arange(num_a) - ind_b = num_a + np.arange(b.shape[1]) + ind_b = num_a + np.arange(b.shape[-1]) # Count number of occurences in combine array, and in a only num_occ_a_and_b = count[k[ind_a]] diff --git a/tests/test_setmembership.py b/tests/test_setmembership.py index e0eb646822..6e8fec56ed 100644 --- a/tests/test_setmembership.py +++ b/tests/test_setmembership.py @@ -111,6 +111,29 @@ def test_ismember_rows_double_occurence_a_and_b(self): assert np.allclose(ma, ma_known) assert np.allclose(ia, ia_known) + def test_ismember_rows_1d(self): + a = np.array([0, 2, 1, 3, 0]) + b = np.array([2, 4, 3]) + + ma, ia = setmembership.ismember_rows(a, b) + + ma_known = np.array([0, 1, 0, 1, 0], dtype=bool) + ia_known = np.array([0, 2]) + + assert np.allclose(ma, ma_known) + assert np.allclose(ia, ia_known) + + def test_ismember_rows_1d(self): + a = np.array([0, 2, 1, 13, 0]) + b = np.array([2, 4, 13, 0]) + + ma, ia = setmembership.ismember_rows(a, b) + + ma_known = np.array([1, 1, 0, 1, 1], dtype=bool) + ia_known = np.array([3, 0, 2, 3]) + + assert np.allclose(ma, ma_known) + assert np.allclose(ia, ia_known) if __name__ == '__main__': unittest.main() From ce1fad07c87c00701364f0801a8cb4ccb5e67222 Mon Sep 17 00:00:00 2001 From: Runar Date: Fri, 24 Mar 2017 11:47:16 +0100 Subject: [PATCH 36/42] Fixed bug in coloring scheme in the case no nodes are connected --- graph.py | 23 ++++++++++++----------- half_space.py | 1 + setmembership.py | 27 ++++++++++++++------------- 3 files changed, 27 insertions(+), 24 deletions(-) diff --git a/graph.py b/graph.py index d2adcc72ed..89fb34b4b3 100644 --- a/graph.py +++ b/graph.py @@ -1,6 +1,7 @@ import numpy as np from scipy import sparse as sps + class Graph: """ A graph class. @@ -22,26 +23,27 @@ class Graph: same colors and nodes in different regions are given different colors. """ - def __init__(self,node_connections): + + def __init__(self, node_connections): self.node_connections = node_connections self.regions = 0 - self.color = np.array([np.NaN]*node_connections.shape[1]) + self.color = np.array([np.NaN] * node_connections.shape[1]) def color_nodes(self): """ Color the nodes in each region """ color = 0 - while self.regions 1: hit = hit[0] ind_of_a_in_b = np.append(ind_of_a_in_b, hit) @@ -149,6 +149,7 @@ def ismember_rows(a, b, sort=True): #--------------------------------------------------------- + def unique_columns_tol(mat, tol=1e-8, exponent=2): """ Remove duplicates from a point set, for a given distance traveling. @@ -164,13 +165,13 @@ def unique_columns_tol(mat, tol=1e-8, exponent=2): the points (due to rounding errors). Defaults to 1e-8. exponent (double, optional): Exponnet in norm used in distance calculation. Defaults to 2. - + Returns: np.ndarray: Unique columns. new_2_old: Index of which points that are preserved old_2_new: Index of the representation of old points in the reduced list. - + Example (won't work as doctest): >>> p_un, n2o, o2n = unique_columns(np.array([[1, 0, 1], [1, 0, 1]])) >>> p_un @@ -188,25 +189,26 @@ def dist(p, pset): pt = p.reshape((-1, 1)) else: pt = p - + return np.power(np.sum(np.power(np.abs(pt - pset), exponent), - axis=0), 1/exponent) + axis=0), 1 / exponent) (nd, l) = mat.shape - + # By default, no columns are kept keep = np.zeros(l, dtype=np.bool) - + # We will however keep the first point keep[0] = True keep_counter = 1 # Map from old points to the unique subspace. Defaults to map to itself. old_2_new = np.arange(l) - + # Loop over all points, check if it is already represented in the kept list for i in range(1, l): - proximate = np.argwhere(dist(mat[:, i], mat[:, keep]) < tol * np.sqrt(nd)) + proximate = np.argwhere( + dist(mat[:, i], mat[:, keep]) < tol * np.sqrt(nd)) if proximate.size > 0: # We will not keep this point @@ -218,6 +220,5 @@ def dist(p, pset): keep_counter += 1 # Finally find which elements we kept new_2_old = np.argwhere(keep).ravel() - - return mat[:, keep], new_2_old, old_2_new + return mat[:, keep], new_2_old, old_2_new From 156a019477670b36d5906bd79ac322f1e7d4aec5 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Fri, 24 Mar 2017 19:19:24 +0100 Subject: [PATCH 37/42] Resolves issue #3 regarding ismember_rows for different types of int. --- setmembership.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/setmembership.py b/setmembership.py index 31ffcf3276..44ca2cbc64 100644 --- a/setmembership.py +++ b/setmembership.py @@ -149,8 +149,11 @@ def ismember_rows(a, b, sort=True, simple_version=False): return_counts=True) else: # Represent the arrays as voids to facilitate quick comparison - voida = _asvoid(sa.transpose()) - voidb = _asvoid(sb.transpose()) + # Take void type of int64s, or else spurious error messages may + # arise. We do this after the transpose (which copies sa and sb) to + # avoid disturbing the original fields. + voida = _asvoid(sa.transpose().astype('int64')) + voidb = _asvoid(sb.transpose().astype('int64')) # Use unique to count the number of occurences in a unq, j, k, count = np.unique(np.vstack((voida, voidb)), From 9b982acb1fcc702473f7e472c8bf952a8800417c Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Sun, 2 Apr 2017 11:35:48 +0200 Subject: [PATCH 38/42] Bugfix in unique_columns_tol in setmembership, related to empty matrices. --- setmembership.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/setmembership.py b/setmembership.py index 44ca2cbc64..111a7e0e37 100644 --- a/setmembership.py +++ b/setmembership.py @@ -223,6 +223,10 @@ def unique_columns_tol(mat, tol=1e-8, exponent=2): """ + # Special treatment of the case with an empty array + if mat.shape[1] == 0: + return mat + def dist(p, pset): " Helper function to compute distance " if p.ndim == 1: From 07dafbcf699cd55c3423ce419426ffe3a63467b2 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Wed, 5 Apr 2017 13:07:13 +0200 Subject: [PATCH 39/42] Unique_columns uses function from numpy 1.13 for integers. Speedup. --- setmembership.py | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/setmembership.py b/setmembership.py index a82e120b9e..f3debd3f76 100644 --- a/setmembership.py +++ b/setmembership.py @@ -7,6 +7,7 @@ import numpy as np +import utils.unique as np_unique def unique_rows(data): """ @@ -226,7 +227,15 @@ def unique_columns_tol(mat, tol=1e-8, exponent=2): # Special treatment of the case with an empty array if mat.shape[1] == 0: - return mat + return mat + + # If the matrix is integers, and the tolerance less than 1/2, we can use + # the new unique function that ships with numpy 1.13. + if issubclass(mat.dtype.type, np.integer) and tol < 0.5: + un_ar, new_2_old, old_2_new \ + = np_unique.unique_np1130(mat, return_index=True, + return_inverse=True, axis=1) + return un_ar, new_2_old, old_2_new def dist(p, pset): " Helper function to compute distance " From edc06637c21f73379387940758f713294dec9b06 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Wed, 5 Apr 2017 13:09:14 +0200 Subject: [PATCH 40/42] Update unit test for unique_columns - returned arrays for integers are now sorted. Necessary because integer input is now delegated to a numpy function. --- tests/test_setmembership.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/tests/test_setmembership.py b/tests/test_setmembership.py index 6e8fec56ed..d63d2045d7 100644 --- a/tests/test_setmembership.py +++ b/tests/test_setmembership.py @@ -148,7 +148,7 @@ def test_no_common_points(self): assert np.allclose(p, p_unique) assert np.alltrue(old_2_new == np.arange(3)) assert np.alltrue(old_2_new == new_2_old) - + def test_remove_one_point(self): p = np.ones((2, 2)) p_unique, new_2_old, old_2_new = setmembership.unique_columns_tol(p) @@ -156,14 +156,14 @@ def test_remove_one_point(self): assert np.allclose(p, np.ones((2, 1))) assert np.alltrue(old_2_new == np.zeros(2)) assert np.alltrue(new_2_old == np.zeros(1)) - + def test_remove_one_of_tree(self): p = np.array([[1, 1, 0], [1, 1, 0]]) p_unique, new_2_old, old_2_new = setmembership.unique_columns_tol(p) - assert np.allclose(p_unique, np.array([[1, 0], [1, 0]])) - assert np.alltrue(old_2_new == np.array([0, 0, 1])) - assert np.alltrue(new_2_old == np.array([0, 2])) + assert np.allclose(p_unique, np.array([[0, 1], [0, 1]])) + assert np.alltrue(old_2_new == np.array([1, 1, 0])) + assert np.alltrue(new_2_old == np.array([2, 0])) if __name__ == '__main__': unittest.main() From a7238c1041f676a40b4b630f087d16cb254b0d2b Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Fri, 21 Apr 2017 20:16:59 +0200 Subject: [PATCH 41/42] Improved comments in utility functions. --- matrix_compression.py | 16 ++++++++------ mcolon.py | 50 +++++++++++++++++++++++++++++++------------ 2 files changed, 45 insertions(+), 21 deletions(-) diff --git a/matrix_compression.py b/matrix_compression.py index 1477fa3862..ff28641a2c 100644 --- a/matrix_compression.py +++ b/matrix_compression.py @@ -1,12 +1,15 @@ -# -*- coding: utf-8 -*- -""" -Created on Tue Mar 1 08:36:47 2016 +""" Functions for compressing matrices to compact format, and recover them. -@author: keile +Acknowledgements: + The functions are a python translation of the corresponding matlab + functions found in the Matlab Reservoir Simulation Toolbox (MRST) developed + by SINTEF ICT, see www.sintef.no/projectweb/mrst/ . + """ import numpy as np + def rldecode(A,n): """ Decode compressed information. @@ -42,7 +45,6 @@ def rlencode(A): return A[::, i], num - if __name__ == '__main__': - A = np.array([[1, 2, 2, 3],[2, 2, 2, 3]]) - rlencode(A) \ No newline at end of file + import doctest + doctest.testmod() diff --git a/mcolon.py b/mcolon.py index 1eef9cb949..9ef700d632 100644 --- a/mcolon.py +++ b/mcolon.py @@ -1,24 +1,41 @@ -# -*- coding: utf-8 -*- -""" -Created on Fri Mar 4 21:02:56 2016 +""" Efficient numpy.arange for arrays of start and end indices. + +Acknowledgements: + The functions are a python translation of the corresponding matlab + functions found in the Matlab Reservoir Simulation Toolbox (MRST) developed + by SINTEF ICT, see www.sintef.no/projectweb/mrst/ . -@author: keile """ import numpy as np def mcolon(lo, hi): - """ Get set of expanded indices - - >>> mcolon(np.array([0, 0, 0]), np.array([1, 3, 2])) - array([0, 1, 0, 1, 2, 3, 0, 1, 2], dtype=int64) - - >>> mcolon(np.array([0, 1]), np.array([2])) - array([0, 1, 2, 1, 2], dtype=int64) - - >>> mcolon(np.array([0, 1, 1, 1]), np.array([0, 3, 3, 3])) - array([0, 1, 2, 3, 1, 2, 3, 1, 2, 3], dtype=int64) + """ Expansion of np.arange(a, b) for arrays a and b. + + The code is equivalent to the following (less efficient) loop: + arr = np.empty(0) + for l, h in zip(lo, hi): + arr = np.hstack((arr, np.arange(l, h+1, 1))) + + Parameters: + lo (np.ndarray, int): Lower bounds of the arrays to be created. + hi (np.ndarray, int): Upper bounds of the arrays to be created. The + elements in hi will be included in the resulting array. + + lo and hi should either have 1 or n elements. If their size are both + larger than one, they should have the same length. + + Examples: + >>> mcolon(np.array([0, 0, 0]), np.array([1, 3, 2])) + array([0, 1, 0, 1, 2, 3, 0, 1, 2], dtype=int64) + + >>> mcolon(np.array([0, 1]), np.array([2])) + array([0, 1, 2, 1, 2], dtype=int64) + + >>> mcolon(np.array([0, 1, 1, 1]), np.array([0, 3, 3, 3])) + array([0, 1, 2, 3, 1, 2, 3, 1, 2, 3], dtype=int64) + """ if lo.size == 1: lo = lo * np.ones(hi.size, dtype='int64') @@ -41,3 +58,8 @@ def mcolon(lo, hi): x[0] = lo[0] x[np.cumsum(d[0:-1]).astype('int64')] = lo[1:] - hi[0:-1] return np.cumsum(x) + + +if __name__ == '__main__': + import doctest + doctest.testmod() From 67191a24523b557d6e8f0ef2f51fd4cd2d334350 Mon Sep 17 00:00:00 2001 From: Eirik Keilegavlen Date: Mon, 24 Apr 2017 12:52:31 +0200 Subject: [PATCH 42/42] Move utils repo to new location --- __init__.py => src/porepy/utils/__init__.py | 0 accumarray.py => src/porepy/utils/accumarray.py | 0 {errors => src/porepy/utils/errors}/error.py | 0 graph.py => src/porepy/utils/graph.py | 0 half_space.py => src/porepy/utils/half_space.py | 0 matrix_compression.py => src/porepy/utils/matrix_compression.py | 0 mcolon.py => src/porepy/utils/mcolon.py | 0 permutations.py => src/porepy/utils/permutations.py | 0 setmembership.py => src/porepy/utils/setmembership.py | 0 {tests => src/porepy/utils/tests}/__init__.py | 0 {tests => src/porepy/utils/tests}/test_graph.py | 0 {tests => src/porepy/utils/tests}/test_half_space.py | 0 {tests => src/porepy/utils/tests}/test_mcolon.py | 0 {tests => src/porepy/utils/tests}/test_permutations.py | 0 {tests => src/porepy/utils/tests}/test_setmembership.py | 0 unique.py => src/porepy/utils/unique.py | 0 16 files changed, 0 insertions(+), 0 deletions(-) rename __init__.py => src/porepy/utils/__init__.py (100%) rename accumarray.py => src/porepy/utils/accumarray.py (100%) rename {errors => src/porepy/utils/errors}/error.py (100%) rename graph.py => src/porepy/utils/graph.py (100%) rename half_space.py => src/porepy/utils/half_space.py (100%) rename matrix_compression.py => src/porepy/utils/matrix_compression.py (100%) rename mcolon.py => src/porepy/utils/mcolon.py (100%) rename permutations.py => src/porepy/utils/permutations.py (100%) rename setmembership.py => src/porepy/utils/setmembership.py (100%) rename {tests => src/porepy/utils/tests}/__init__.py (100%) rename {tests => src/porepy/utils/tests}/test_graph.py (100%) rename {tests => src/porepy/utils/tests}/test_half_space.py (100%) rename {tests => src/porepy/utils/tests}/test_mcolon.py (100%) rename {tests => src/porepy/utils/tests}/test_permutations.py (100%) rename {tests => src/porepy/utils/tests}/test_setmembership.py (100%) rename unique.py => src/porepy/utils/unique.py (100%) diff --git a/__init__.py b/src/porepy/utils/__init__.py similarity index 100% rename from __init__.py rename to src/porepy/utils/__init__.py diff --git a/accumarray.py b/src/porepy/utils/accumarray.py similarity index 100% rename from accumarray.py rename to src/porepy/utils/accumarray.py diff --git a/errors/error.py b/src/porepy/utils/errors/error.py similarity index 100% rename from errors/error.py rename to src/porepy/utils/errors/error.py diff --git a/graph.py b/src/porepy/utils/graph.py similarity index 100% rename from graph.py rename to src/porepy/utils/graph.py diff --git a/half_space.py b/src/porepy/utils/half_space.py similarity index 100% rename from half_space.py rename to src/porepy/utils/half_space.py diff --git a/matrix_compression.py b/src/porepy/utils/matrix_compression.py similarity index 100% rename from matrix_compression.py rename to src/porepy/utils/matrix_compression.py diff --git a/mcolon.py b/src/porepy/utils/mcolon.py similarity index 100% rename from mcolon.py rename to src/porepy/utils/mcolon.py diff --git a/permutations.py b/src/porepy/utils/permutations.py similarity index 100% rename from permutations.py rename to src/porepy/utils/permutations.py diff --git a/setmembership.py b/src/porepy/utils/setmembership.py similarity index 100% rename from setmembership.py rename to src/porepy/utils/setmembership.py diff --git a/tests/__init__.py b/src/porepy/utils/tests/__init__.py similarity index 100% rename from tests/__init__.py rename to src/porepy/utils/tests/__init__.py diff --git a/tests/test_graph.py b/src/porepy/utils/tests/test_graph.py similarity index 100% rename from tests/test_graph.py rename to src/porepy/utils/tests/test_graph.py diff --git a/tests/test_half_space.py b/src/porepy/utils/tests/test_half_space.py similarity index 100% rename from tests/test_half_space.py rename to src/porepy/utils/tests/test_half_space.py diff --git a/tests/test_mcolon.py b/src/porepy/utils/tests/test_mcolon.py similarity index 100% rename from tests/test_mcolon.py rename to src/porepy/utils/tests/test_mcolon.py diff --git a/tests/test_permutations.py b/src/porepy/utils/tests/test_permutations.py similarity index 100% rename from tests/test_permutations.py rename to src/porepy/utils/tests/test_permutations.py diff --git a/tests/test_setmembership.py b/src/porepy/utils/tests/test_setmembership.py similarity index 100% rename from tests/test_setmembership.py rename to src/porepy/utils/tests/test_setmembership.py diff --git a/unique.py b/src/porepy/utils/unique.py similarity index 100% rename from unique.py rename to src/porepy/utils/unique.py