From 3c7839891492c639b597c8b18860ca849b2988b4 Mon Sep 17 00:00:00 2001 From: Gabriel Girard Date: Fri, 20 Jan 2023 11:25:46 -0500 Subject: [PATCH 1/8] NF - added cython utility fcts --- dipy/utils/fast_numpy.pxd | 25 +++++++++++ dipy/utils/fast_numpy.pyx | 90 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 115 insertions(+) diff --git a/dipy/utils/fast_numpy.pxd b/dipy/utils/fast_numpy.pxd index 6d4943b0fe..e8bd7d4b91 100644 --- a/dipy/utils/fast_numpy.pxd +++ b/dipy/utils/fast_numpy.pxd @@ -4,6 +4,13 @@ cimport numpy as np +from libc.stdlib cimport rand +from libc.math cimport sqrt, fabs, M_PI, pow, sin, cos + + +cdef extern from "limits.h": + int INT_MAX + # Replaces a numpy.searchsorted(arr, number, 'right') cdef int where_to_insert( np.float_t* arr, @@ -22,3 +29,21 @@ cdef void copy_point( cdef void scalar_muliplication_point( double * a, double scalar) nogil + +cdef double random() + +cdef double norm( + double[:] v) + +cdef double dot( + double[:] v1, + double[:] v2) + +cdef void normalize( + double[:] v) + +cdef void cross( + double[:] out, + double[:] v1, + double[:] v2) + \ No newline at end of file diff --git a/dipy/utils/fast_numpy.pyx b/dipy/utils/fast_numpy.pyx index f0a577f949..1b9fd10fca 100644 --- a/dipy/utils/fast_numpy.pyx +++ b/dipy/utils/fast_numpy.pyx @@ -45,3 +45,93 @@ cdef void scalar_muliplication_point( int i = 0 for i in range(3): a[i] *= scalar + + +cdef double random(): + """Sample a random number between 0 and 1. + + Returns + ------- + _ : double + random number + """ + return rand() / float(INT_MAX) + + +cdef double norm(double[:] v): + """Vector norm. + + Parameters + ---------- + v : double[3] + input vector. + + Returns + ------- + _ : double + norm of the vector. + """ + return sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]) + + +cdef double dot(double[:] v1, double[:] v2): + """Dot product. + + Parameters + ---------- + v1 : double[3] + input vector 1. + v2 : double[3] + input vector 2. + + Returns + ------- + _ : double + dot product of input vectors. + """ + return v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2] + + +cdef void normalize(double[:] v): + """Vector normalization. + + Parameters + ---------- + v : double[3] + input vector + + Notes + ----- + Overwrites the first argument. + + """ + cdef double scale = 1.0 / norm(v) + v[0] = v[0] * scale + v[1] = v[1] * scale + v[2] = v[2] * scale + + + + + +cdef void cross(double[:] out, double[:] v1, double[:] v2): + """Cross product. + + Parameters + ---------- + out : double[3] + output. + v1 : double[3] + input vector 1. + v2 : double[3] + input vector 2. + + Notes + ----- + Overwrites the first argument. + + + """ + out[0] = v1[1] * v2[2] - v1[2] * v2[1] + out[1] = v1[2] * v2[0] - v1[0] * v2[2] + out[2] = v1[0] * v2[1] - v1[1] * v2[0] \ No newline at end of file From 83e4a34c338bfdeeab10a325c1c134219ea6e84e Mon Sep 17 00:00:00 2001 From: Gabriel Girard Date: Fri, 20 Jan 2023 11:29:40 -0500 Subject: [PATCH 2/8] RF - removed unused impots --- dipy/utils/fast_numpy.pxd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dipy/utils/fast_numpy.pxd b/dipy/utils/fast_numpy.pxd index e8bd7d4b91..3e0114db01 100644 --- a/dipy/utils/fast_numpy.pxd +++ b/dipy/utils/fast_numpy.pxd @@ -5,7 +5,7 @@ cimport numpy as np from libc.stdlib cimport rand -from libc.math cimport sqrt, fabs, M_PI, pow, sin, cos +from libc.math cimport sqrt cdef extern from "limits.h": From 05dd7a39c904818c57071adbdec058472c545e20 Mon Sep 17 00:00:00 2001 From: Gabriel Girard Date: Fri, 20 Jan 2023 11:56:14 -0500 Subject: [PATCH 3/8] DOC/PEP8 --- dipy/utils/fast_numpy.pyx | 19 ++++++++----------- 1 file changed, 8 insertions(+), 11 deletions(-) diff --git a/dipy/utils/fast_numpy.pyx b/dipy/utils/fast_numpy.pyx index 1b9fd10fca..61eed85988 100644 --- a/dipy/utils/fast_numpy.pyx +++ b/dipy/utils/fast_numpy.pyx @@ -59,7 +59,7 @@ cdef double random(): cdef double norm(double[:] v): - """Vector norm. + """Compute the vector norm. Parameters ---------- @@ -69,13 +69,13 @@ cdef double norm(double[:] v): Returns ------- _ : double - norm of the vector. + norm of the input vector. """ return sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]) cdef double dot(double[:] v1, double[:] v2): - """Dot product. + """Compute vectors dot product. Parameters ---------- @@ -93,7 +93,7 @@ cdef double dot(double[:] v1, double[:] v2): cdef void normalize(double[:] v): - """Vector normalization. + """Normalize the vector. Parameters ---------- @@ -111,16 +111,13 @@ cdef void normalize(double[:] v): v[2] = v[2] * scale - - - cdef void cross(double[:] out, double[:] v1, double[:] v2): - """Cross product. + """Compute vectors cross product. Parameters ---------- out : double[3] - output. + output vector. v1 : double[3] input vector 1. v2 : double[3] @@ -129,9 +126,9 @@ cdef void cross(double[:] out, double[:] v1, double[:] v2): Notes ----- Overwrites the first argument. - """ out[0] = v1[1] * v2[2] - v1[2] * v2[1] out[1] = v1[2] * v2[0] - v1[0] * v2[2] - out[2] = v1[0] * v2[1] - v1[1] * v2[0] \ No newline at end of file + out[2] = v1[0] * v2[1] - v1[1] * v2[0] + \ No newline at end of file From 64e2f394d4948ac2bd9ca5d25487dd8b669ee912 Mon Sep 17 00:00:00 2001 From: Gabriel Girard Date: Fri, 20 Jan 2023 15:13:30 -0500 Subject: [PATCH 4/8] RF - change cdef to cpdef for tests --- dipy/utils/fast_numpy.pxd | 18 +++++++++--------- dipy/utils/fast_numpy.pyx | 10 +++++----- 2 files changed, 14 insertions(+), 14 deletions(-) diff --git a/dipy/utils/fast_numpy.pxd b/dipy/utils/fast_numpy.pxd index 3e0114db01..914af00db7 100644 --- a/dipy/utils/fast_numpy.pxd +++ b/dipy/utils/fast_numpy.pxd @@ -30,20 +30,20 @@ cdef void scalar_muliplication_point( double * a, double scalar) nogil -cdef double random() +cpdef double random() nogil -cdef double norm( - double[:] v) +cpdef double norm( + double[:] v) nogil -cdef double dot( +cpdef double dot( double[:] v1, - double[:] v2) + double[:] v2) nogil -cdef void normalize( - double[:] v) +cpdef void normalize( + double[:] v) nogil -cdef void cross( +cpdef void cross( double[:] out, double[:] v1, - double[:] v2) + double[:] v2) nogil \ No newline at end of file diff --git a/dipy/utils/fast_numpy.pyx b/dipy/utils/fast_numpy.pyx index 61eed85988..40ec4c01d3 100644 --- a/dipy/utils/fast_numpy.pyx +++ b/dipy/utils/fast_numpy.pyx @@ -47,7 +47,7 @@ cdef void scalar_muliplication_point( a[i] *= scalar -cdef double random(): +cpdef double random() nogil: """Sample a random number between 0 and 1. Returns @@ -58,7 +58,7 @@ cdef double random(): return rand() / float(INT_MAX) -cdef double norm(double[:] v): +cpdef double norm(double[:] v) nogil: """Compute the vector norm. Parameters @@ -74,7 +74,7 @@ cdef double norm(double[:] v): return sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]) -cdef double dot(double[:] v1, double[:] v2): +cpdef double dot(double[:] v1, double[:] v2) nogil: """Compute vectors dot product. Parameters @@ -92,7 +92,7 @@ cdef double dot(double[:] v1, double[:] v2): return v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2] -cdef void normalize(double[:] v): +cpdef void normalize(double[:] v) nogil: """Normalize the vector. Parameters @@ -111,7 +111,7 @@ cdef void normalize(double[:] v): v[2] = v[2] * scale -cdef void cross(double[:] out, double[:] v1, double[:] v2): +cpdef void cross(double[:] out, double[:] v1, double[:] v2) nogil: """Compute vectors cross product. Parameters From 956a2bfa8d2e096ee10350914b71358310070794 Mon Sep 17 00:00:00 2001 From: Gabriel Girard Date: Fri, 20 Jan 2023 15:14:00 -0500 Subject: [PATCH 5/8] TST - fast numpy fcts --- dipy/utils/tests/test_fast_numpy.py | 75 +++++++++++++++++++++++++++++ 1 file changed, 75 insertions(+) create mode 100644 dipy/utils/tests/test_fast_numpy.py diff --git a/dipy/utils/tests/test_fast_numpy.py b/dipy/utils/tests/test_fast_numpy.py new file mode 100644 index 0000000000..78324b9755 --- /dev/null +++ b/dipy/utils/tests/test_fast_numpy.py @@ -0,0 +1,75 @@ +import timeit + +import numpy as np + +from numpy.testing import (assert_array_equal, assert_array_almost_equal, + assert_almost_equal, assert_equal, assert_) +from dipy.utils.fast_numpy import random, norm, normalize, dot, cross + + +def test_random(): + # Test that random numbers are between 0 and 1. + for _ in range(10): + vec = random() + assert_(vec <= 1 and vec >= 0) + + +def test_norm(): + # Test that the norm equal numpy norm. + for _ in range(10): + vec = np.random.random(3) + assert_equal(norm(vec), np.linalg.norm(vec)) + + +def test_normalize(): + # Test that the normalize vector as a norm of 1. + for _ in range(10): + vec = np.random.random(3) + normalize(vec) + assert_almost_equal(np.linalg.norm(vec), 1) + +def test_dot(): + # Test that dot is faster and equal to numpy.dot. + for _ in range(10): + vec1 = np.random.random(3) + vec2 = np.random.random(3) + assert_equal(dot(vec1, vec2), np.dot(vec1, vec2)) + + vec1 = np.random.random(3) + vec2 = np.random.random(3) + + def __dot(): + dot(vec1, vec2) + + def __npdot(): + np.dot(vec1, vec2) + + number=100000 + time_dot = timeit.timeit(__dot, number=number) + time_npdot = timeit.timeit(__npdot, number=number) + assert_(time_dot < time_npdot) + + +def test_cross(): + # Test that cross is faster and equal to numpy.cross. + out = np.zeros(3) + for _ in range(10): + vec1 = np.random.random(3) + vec2 = np.random.random(3) + cross(out, vec1, vec2) + assert_array_equal(out, np.cross(vec1, vec2)) + + vec1 = np.random.random(3) + vec2 = np.random.random(3) + + def __cross(): + cross(out, vec1, vec2) + + def __npcross(): + np.cross(vec1, vec2) + + number=10000 + time_cross = timeit.timeit(__cross, number=number) + time_npcross = timeit.timeit(__npcross, number=number) + assert_(time_cross < time_npcross) + print(time_cross, time_npcross) From c7a8a31b1142c82a0281b3c0b2f08c87f5c82673 Mon Sep 17 00:00:00 2001 From: Gabriel Girard Date: Fri, 20 Jan 2023 15:18:17 -0500 Subject: [PATCH 6/8] RF - removed prints --- dipy/utils/tests/test_fast_numpy.py | 1 - 1 file changed, 1 deletion(-) diff --git a/dipy/utils/tests/test_fast_numpy.py b/dipy/utils/tests/test_fast_numpy.py index 78324b9755..90cfa546c0 100644 --- a/dipy/utils/tests/test_fast_numpy.py +++ b/dipy/utils/tests/test_fast_numpy.py @@ -72,4 +72,3 @@ def __npcross(): time_cross = timeit.timeit(__cross, number=number) time_npcross = timeit.timeit(__npcross, number=number) assert_(time_cross < time_npcross) - print(time_cross, time_npcross) From 7ef818057acbd6ddf9a8e17b8642fa77865d577a Mon Sep 17 00:00:00 2001 From: Gabriel Girard Date: Thu, 26 Jan 2023 13:16:23 -0500 Subject: [PATCH 7/8] TST - changed assert_equal to assert_almost_equal --- dipy/utils/tests/test_fast_numpy.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/dipy/utils/tests/test_fast_numpy.py b/dipy/utils/tests/test_fast_numpy.py index 90cfa546c0..74bf40f32a 100644 --- a/dipy/utils/tests/test_fast_numpy.py +++ b/dipy/utils/tests/test_fast_numpy.py @@ -2,8 +2,7 @@ import numpy as np -from numpy.testing import (assert_array_equal, assert_array_almost_equal, - assert_almost_equal, assert_equal, assert_) +from numpy.testing import assert_almost_equal, assert_ from dipy.utils.fast_numpy import random, norm, normalize, dot, cross @@ -18,7 +17,7 @@ def test_norm(): # Test that the norm equal numpy norm. for _ in range(10): vec = np.random.random(3) - assert_equal(norm(vec), np.linalg.norm(vec)) + assert_almost_equal(norm(vec), np.linalg.norm(vec)) def test_normalize(): @@ -33,7 +32,7 @@ def test_dot(): for _ in range(10): vec1 = np.random.random(3) vec2 = np.random.random(3) - assert_equal(dot(vec1, vec2), np.dot(vec1, vec2)) + assert_almost_equal(dot(vec1, vec2), np.dot(vec1, vec2)) vec1 = np.random.random(3) vec2 = np.random.random(3) From 9393354aa2b925707a90960cde72fd84d8623a51 Mon Sep 17 00:00:00 2001 From: Gabriel Girard Date: Thu, 26 Jan 2023 15:21:20 -0500 Subject: [PATCH 8/8] TST/BF - add missing import --- dipy/utils/tests/test_fast_numpy.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/dipy/utils/tests/test_fast_numpy.py b/dipy/utils/tests/test_fast_numpy.py index 74bf40f32a..1ab04e28c5 100644 --- a/dipy/utils/tests/test_fast_numpy.py +++ b/dipy/utils/tests/test_fast_numpy.py @@ -1,8 +1,7 @@ import timeit import numpy as np - -from numpy.testing import assert_almost_equal, assert_ +from numpy.testing import assert_, assert_almost_equal, assert_array_equal from dipy.utils.fast_numpy import random, norm, normalize, dot, cross