Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

ENH: added vector operators: divergence, curl and laplacian #6727

Closed
wants to merge 8 commits into from
Closed

ENH: added vector operators: divergence, curl and laplacian #6727

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
517a053
ENH: added vector operators: divergence, curl and laplacian
Nov 25, 2015
57d1147
ENH: added vector operators: divergence, curl and laplacian
Nov 26, 2015
dfae939
ENH: added vector operators: divergence, curl and laplacian
Nov 29, 2015
e673560
ENH: added vector operators: divergence, curl and laplacian
Nov 30, 2015
8ac65a7
ENH: added vector operators: divergence, curl and laplacian
Dec 2, 2015
87c19b8
ENH: added vector operators: divergence, curl and laplacian
Dec 2, 2015
e71dbcd
ENH: added vector operators: divergence, curl and laplacian
Dec 20, 2015
eeada1f
ENH: added vector operators: divergence, curl and laplacian
Dec 20, 2015
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
20 changes: 20 additions & 0 deletions doc/release/1.11.0-notes.rst 100644 → 100755
View file
Open in desktop
Expand Up @@ -89,6 +89,26 @@ The function now internally calls the generic ``npy_amergesort``
when the type does not implement a merge-sort kind of ``argsort``
method.

New function *divergence*
~~~~~~~~~~~~~~~~~~~~~~~~~
*divergence* makes use of the gradient function to calculate the divergence
of a list of arrays (vector field) as the sum of the derivatives of each
component with respect to the corresponding axis.

New function *curl*
~~~~~~~~~~~~~~~~~~~
*curl* makes use of gradient to calculate the curl of lists of 3 arrays
(3-D vector fields).

New function *laplace*
~~~~~~~~~~~~~~~~~~~~~~
*laplace* calculates the laplacian of either an array (scalar laplacian) or
a list of arrays (vector field). In the first case, a composition of the
functions gradient and divergence is used to return an array. In the second
case, the same operation is applied to every component (array) independently,
and a list of these is returned.


Changes
=======

Expand Down
211 changes: 210 additions & 1 deletion numpy/lib/function_base.py 100644 → 100755
View file
Open in desktop
Expand Up @@ -41,7 +41,8 @@
'histogram', 'histogramdd', 'bincount', 'digitize', 'cov', 'corrcoef',
'msort', 'median', 'sinc', 'hamming', 'hanning', 'bartlett',
'blackman', 'kaiser', 'trapz', 'i0', 'add_newdoc', 'add_docstring',
'meshgrid', 'delete', 'insert', 'append', 'interp', 'add_newdoc_ufunc'
'meshgrid', 'delete', 'insert', 'append', 'interp', 'add_newdoc_ufunc',
'_inplace_sum', 'divergence', 'curl', 'laplace'
]


Expand Down Expand Up @@ -1320,6 +1321,214 @@ def gradient(f, *varargs, **kwargs):
else:
return outvals

def _inplace_sum(items):
"""
Returns the sum of the elements of a list.
Used in the follwing functions:
divergence,

.. versionadded:: 1.11.0

"""
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Blank line before """.

it = iter(items)
total = next(it)
for item in it:
total += item
return total
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this can be written np.sum(arr, out=arr[0]). Seems like a premature optimization to me though



def divergence(v,*varargs):
"""
Return the divergence of an N-dimensional vector field of N
components, each of dimension N.

The divergence is computed using second order accurate central
differences in the interior and either first differences or second
order accurate one-sides (forward or backwards) differences at the
boundaries. The returned gradient hence has the same shape as the
input array.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should have .. versionadded:: 1.11.0 here. Grep for examples.

Parameters
----------
v : list of numpy arrays
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Wouldn't it make more sense to make this an array-like? Here, one can choose which dimension to be that of the axes; I think the first one is fine. In that case, one could just do

v = np.asanyarray(v)

and it would even cover the case of having a list of numpy arrays. Another advantage would be that the tests for the sizes to be the same could be omitted, since that would be done implicitly already.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi,
Thanks for your suggestion. By reviewing the size test, I have realized that they are incomplete right now, they only check that every component has the same number of dimensions, but not the same size. A quick example:

a = np.array([[1,2],[3,4]])
b = np.array([[5,6,7],[8,9,10],[11,12,13]])
c=[a,b]
d = np.asanyarray(c)
np.shape(d[0])
(2,2)
np.shape(d[1])
(3,3)
This should fail the tests for size. But the tests right now would pass since:
len(d[0].shape)
2
len(d[1].shape)
2

I'm going to fix this and will implement your suggestion as well.
Thank you

Each of these array is the N-dimensional projection of the vector
field on the corresponding axis.

varargs : scalar or list of scalar, optional
N scalars specifying the sample distances for each dimension,
i.e. `dx`, `dy`, `dz`, ... Default distance: 1.
single scalar specifies sample distance for all dimensions.
if `axis` is given, the number of varargs must equal the number of axes.

.. versionadded:: 1.11.0

Returns
-------
divergence : numpy array
The output corresponds to the divergence of the input vector field.
This means that the output array has the form
dAx/dx + dAy/dy + dAz/dz + ... for an input vector field
A = (Ax, Ay, Az, ...) up to N dimensions.

Example
-------
This example shows how the calculated divergence of the 2-D field
(0.5*x**2, -y*x) (whose divergence should be 0 everywhere) returns a 0 array

>>> import numpy as np
>>> X,Y=np.mgrid[0:2000,0:2000]
>>> a1=0.5*X**2
>>> a2=-Y*X
>>> c=[a1,a2]
>>> d=np.divergence(c)
>>> print d
[[ 0. 0. 0. ..., 0. 0. 0.]
[ 0. 0. 0. ..., 0. 0. 0.]
[ 0. 0. 0. ..., 0. 0. 0.]
...,
[ 0. 0. 0. ..., 0. 0. 0.]
[ 0. 0. 0. ..., 0. 0. 0.]
[ 0. 0. 0. ..., 0. 0. 0.]

"""

N = len(v)
axes = tuple(range(N))
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this a mistake? Given this, why is there a need to normalize the axes below?

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, you are right, I borrowed some code from np.gradient and forgot to remove the normalization below. I will do it now, thanks.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This seems overly long as well. Doesn't the below cover it?

outvals = gradient(v[0], *varargs, axis=0, edge_order=0)
for axis in range(1, len(v)):
    outvals += gradient(v[axis], *varargs, axis=axis, edge_order=0)
return outvals

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

p.s. Is there any reason not just to have edge_order as a keyword argument and pass it on?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, I was thinking along similar lines.

A trick you could use to avoid repeating the first line outside the loop is to make a dummy object to add in-place:

class ZeroArray(object):
    def __iadd__(self, other):
        return other

The first time you use ZeroArray, it's replaced by the object you add in-place:

>>> x = ZeroArray()
>>> x += 10
>>> x
10
>>> type(x)
int

Then you could write:

outvals = ZeroArray()
for axis in range(len(v)):
    outvals += ...

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Neat. Though I fear it fails the Zen of Python...

Actually, I think python's sum is smart enough to do this intrinsically, i.e., one would just do

outvals = sum(gradient(v[axis], *varargs, axis=axis, edge_order=0)
              for axis in range(len(v)))

This definitely beats everything else in readability, but, doing a quick test, is quite a bit slower than either of our approaches (which both are faster than first making a zero-filled array for large sizes)

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Indeed. The problem with sum is that it doesn't do in-place arithmetic on the result, e.g.,

class A(object):
    def __iadd__(self, other):
        print('iadd with %r' % other)
        return self

    def __add__(self, other):
        print('add with %r' % other)
        return self

    def __radd__(self, other):
        print('radd with %r' % other)
        return self

sum([A()] * 4)
# radd with 0
# add with <__main__.A object at 0x107590748>
# add with <__main__.A object at 0x107590748>
# add with <__main__.A object at 0x107590748>

Maybe defining inplace_sum would be cleaner?

def inplace_sum(items, start=0):
    total = start
    for item in items:
        total += item
    return total

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Actually, in this case we want the result to be added in-place on the first object to avoid an unnecessary memory allocation. So really, we need something like this:

def inplace_sum(items):
    it = iter(items)
    total = next(it)
    for item in it:
        total += item
    return total

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, that makes sense -- I can see how for a generic sum you cannot rely on inplace addition doing what you'd like. You solution also makes sense, and is similarly fast to our other options. Still, probably getting a bit off-topic here -- I think for the divergence I would just stick with doing it explicitly as in my original suggestion.

out = _inplace_sum(gradient(v[ax], *varargs, axis=ax, edge_order=0)
for ax in axes)
return out

def curl(v,*varargs):
"""
Return the curl of a 3-D vector field.

The curl is computed using second order accurate central differences
in the interior and either first differences or second order accurate
one-sides (forward or backwards) differences at the boundaries. The
returned gradient hence has the same shape as the input array.

Parameters
----------
v : list of numpy arrays
Each of these array is the N-dimensional projection
of the vector field on the corresponding axis.

varargs : scalar or list of scalar, optional
N scalars specifying the sample distances for each dimension,
i.e. `dx`, `dy`, `dz`, ... Default distance: 1.
single scalar specifies sample distance for all dimensions.
if `axis` is given, the number of varargs must equal the number of axes.

.. versionadded:: 1.11.0

Returns
-------
curl : list of numpy arrays
The output corresponds to the curl of the input 3-D vector field.
This means that the output has the form
(dAz/dy-dAy/dz, dAx/dz-dAz/dx, dAy/dx-dAx/dy) of an input vector
field A = (Ax, Ay, Ax).

Example
-------
The following example shows in matplotlib how the vector
field (0,0,x*y) has the correct value, nonzero only in the
third component of the resulting vector field after applying curl

>>> import matplotlib.pylab as plt
>>> import numpy as np
>>> X,Y,Z=np.mgrid[0:200,0:200,0:200]
>>> a0=X*Y
>>> a1=np.zeros(np.shape(X))
>>> a2=np.zeros(np.shape(X))
>>> a=[a0,a1,a2]
>>> [cx,cy,cz]=np.curl(a)
>>> plt.imshow(cx[:,:,1])
>>> plt.colorbar()
>>> plt.show()
>>> plt.imshow(cy[:,:,1])
>>> plt.colorbar()
>>> plt.show()
>>> plt.imshow(cz[:,:,1])
>>> plt.colorbar()
>>> plt.show()

"""

if len(v) != 3:
raise ValueError("Enter a list of 3 arrays")

outvals = [gradient(v[(ax+2) % 3], *varargs, axis=(ax+1) % 3,
edge_order=0) -
gradient(v[(ax+1) % 3], *varargs, axis=(ax+2) % 3,
edge_order=0)
for ax in range(3)]
return outvals

def laplace(f, *varargs):
"""
Return the laplacian of input.

Computes the laplacian of an N-dimensional numpy array (scalar
field), or of a list of N N-dimensional numpy arrays (vector
field) as the composition of the numpy functions gradient and
divergence. In the case of an input vector field, the result
is the vector laplacian, in which a list of the laplacian of
each component is returned.
For more information on the computations, type help(gradient) or
help(divergence)

Parameters
----------
f : array_like
Input array

.. versionadded:: 1.11.0

Returns
-------
l : array_like
Laplacian of input. If the input is a scalar field f
(single numpy array) the output is (d^2(f)/dx_0^2 +
d^2(f)/dx_1^2 + d^2(f)/dx_2^2 + ...), where d^2(f)/dx_i^2
refers to the second derivative of f with respect to x_i.
If the input is a scalar field the output is a list of
numpy arrays, being each of them the laplacian operator
applied to each of the input vector field components f_i

See Also
--------
gradient, divergence

Example
-------
This example returns an array that grows linearly in the X axis
after applying the laplacian function to the array X**3+Y**2+Z**2:

>>> X,Y,Z=np.mgrid[0:200,0:200,0:200]
>>> f=X**3+Y**2+Z**2
>>> d=np.laplace(f)
>>> plt.imshow(d[:,:,1])
>>> plt.colorbar()
>>> plt.show()

"""
if type(f) == np.ndarray:
l = divergence(gradient(f,*varargs, edge_order=0),*varargs)
elif type(f) == list:
if not all([np.shape(f[0]) == np.shape(f[i]) for i in range(len(f))]):
raise TypeError("All components of the input vector field "
"must be the same shape")
else:
l = []
for i in range(len(f)):
l.append(divergence(gradient(f[i],
*varargs, edge_order=0)),*varargs)

else:
raise TypeError("Please, enter a numpy array or a list of"
" numpy arrays of the same shape")
return l

def diff(a, n=1, axis=-1):
"""
Expand Down
91 changes: 91 additions & 0 deletions numpy/lib/tests/test_function_base.py 100644 → 100755
View file
Open in desktop
Expand Up @@ -627,6 +627,97 @@ def test_specific_axes(self):
assert_raises(TypeError, gradient, x, axis=[1,])


class TestDivergence(TestCase):

def test_basic(self):
v = [[[1, 2], [3, 4]], [[1, 3], [2, 4]]]
x = np.array(v)
dx = np.array([[4., 4.], [4., 4.]])
assert_array_equal(np.divergence(x), dx)
assert_array_equal(np.divergence(v), dx)

def test_badargs(self):
# for 2D array, divergence can take 0, 1, or 2 extra args
x = np.array([[[1, 2], [3, 4]], [[1, 3], [2, 4]]])
assert_raises(SyntaxError, np.divergence, x, np.array([1., 1.]),
np.array([1., 1.]), np.array([1., 1.]))

def test_masked(self):
# Make sure that divergence supports subclasses like masked arrays
x = np.ma.array([[[1, 1], [3, 4]], [[1, 2], [3, 4]]],
mask=[[[False, False], [False, False]],
[[False, False], [True, False]]])
out = np.divergence(x)
assert_equal(type(out), type(x))
# And make sure that the output and input don't have aliased mask
# arrays
assert_(x.mask is not out.mask)

def test_datetime64(self):
# Make sure that divergence() can handle special
# data types like datetime64
x = np.array(
[[['1910-08-16', '1910-08-11'], ['1910-08-10', '1910-08-12']],
[['1910-08-13', '1910-08-17'], ['1910-08-03', '1910-08-05']]],
dtype='datetime64[D]')
divx = np.array([[-2, 5], [-4, 3]], dtype='timedelta64[D]')
assert_array_equal(np.divergence(x), divx)
assert_(divx.dtype == np.dtype('timedelta64[D]'))

def test_timedelta64(self):
# Make sure gradient() can handle special types like timedelta64
x = np.array([[[1,2], [4,5]], [[1,2], [4,5]]],
dtype = 'timedelta64[D]')
divx = np.array([[4, 4], [4, 4]], dtype='timedelta64[D]')

assert_array_equal(np.divergence(x), divx)
assert_(divx.dtype == np.dtype('timedelta64[D]'))

class TestCurl(TestCase):

def test_basic(self):
a = np.array([[[1, 2], [3, 4]], [[5, 6], [7, 8]]])
x = [a, a, a]
v = np.array(x)
curlx = [np.array([[[1., 1.],[1., 1.]], [[1., 1.],
[1., 1.]]]), np.array([[[-3., -3.], [-3., -3.]], [[-3., -3.],
[-3., -3.]]]), np.array([[[2., 2.], [2., 2.]], [[2., 2.], [2., 2.]]])]
assert_array_equal(np.curl(x), curlx)
assert_array_equal(np.curl(v), curlx)

class TestLaplace(TestCase):

def test_basic_scalar(self):
a = np.array([[3, 1, 5], [6, 8, 14], [2, 9, 14]])
x = np.laplace(a)
lapx = np.array([[-1, 0, -3], [-3, -2, -5], [-9, -8, -11]])
assert_array_equal(x, lapx)

def test_basic_vector(self):
p=np.array([[[1,4,5],[3,6,2],[6,5,4]],[[7,6,5],[1,2,3],[6,9,7]],
[[2,3,5],[7,6,5],[9,7,6]]])
q=np.array([[[1,3,7],[2,6,8],[3,5,1]],[[1,1,5],[1,4,3],[7,4,8]],
[[2,7,55],[1,3,5],[4,1,1]]])
r=np.array([[[2,2,2],[1,1,3],[3,3,2]],[[5,5,5],[7,4,2],[6,6,6]],
[[11,5,25],[10,13,5],[4,11,21]]])

a = [[[-12., -10., 3.], [2., -2., -1.], [4., -9., 1.]],
[[0., 6., 6.], [19., 19., 7.], [9., 0., -3.]],
[[-13., -6., 2.], [5., 6., 2.], [1., -7., -2.]]]

b = [[[3., 6., 46.], [-1., -5., -3.], [-13., -12., -28.]],
[[11., 9., 63.], [3., -6., 10.], [6., 2., 0.]],
[[48., 53., 141.], [5., 3., 53.], [0., 3., 35.]]]

c = [[[6., 0., 15.], [2., 11., 4.], [-3., 4., 8.]],
[[0., 0., 24.], [-5., 10., 12.], [-8., 5., 18.]],
[[24., 13., 79.], [-19.,-15., 29.], [-7., -5., 50.]]]

p1 = [p, q, r]
a1 = [a, b, c]

assert_array_equal(np.laplace(p1), a1)

class TestAngle(TestCase):

def test_basic(self):
Expand Down