making a function that can take arguments in various shapes

Question

Q1) Numpy functions can take arguments in different shapes. For instance, np.sum(V) can take either of two below and return outputs with different shapes.

x1= np.array( [1,3] ) #(1)
x2= np.array([[[1,2],[3,4]], [[5,6],[7,8]]]) #(2)

I am making my own function something like below, which adds two values in an 1D vector with the length of two and return the real number.

def foo(V):
    return V[0]+V[1];

However, this foo function can only take one 1D vector and cannot take any other shapes. It can only take x1 above as an argument but not x2. If I want to make my function work with either of two variables above(x1 and x2), or with any other shapes that has arrays with the length of 2 in their last dimension, how should I revise my foo function?

---

---------------------------update------------------------------

My original function was a hardcoded negative gaussian pdf function.

def nGauss(X, mu, cov):
    # multivariate negative gaussian.    
    # mu is a vector and cov is a covariance matrix.

    k = X.shape[0];
    dev = X-mu
    p1 = np.power( np.power(np.pi * 2, k) , -0.5);
    p2 = np.power( np.linalg.det(cov)  , -0.5)
    p3 = np.exp( -0.5 * np.dot( np.dot(dev.transpose(), np.linalg.inv(cov)), dev));

    return -1.0 * p1 * p2 * p3;

Now his function can return only one pdf value. For example, it can only take arguments like np.array([1,2]), but cannot take arguments X like np.array([[[1,2], [5,6]], [[7,8],[9,0]]]). Here my question was how to make my gaussian function takes arguments of arbitrary shapes and return the pdf value of each point maintaining the same structure except the last dimension, such as nGauss(np.array( [1,2] ), mu, cov) returns [ 0.000023 ], and nGauss(np.array([[[1,2], [5,6]], [[7,8],[9,0]]]), mu, cov) returns [[ 0.000023, 0000014], [0.000012, 0.000042]].

I notice that scipy function 'multivariate_normal.pdf' can do this.

---

Q2) I am also having a difficulty in understanding np's basic array.

t1=np.array([[1,2,3], [4,5,6]])
t2=np.array([1,2,3])
t3=np.array([[1,2,3], [4,5,6],5])

The shape of t1 is (2,3), and it seems legitimate in terms of matrix perspective; 2 rows and 3 columns. However, the shape of t2 is (3,), which I think has to be (3). What's the meaning of the empty space after "3,"? Also, the shape of t3 is (3,). In this case, is the meaning of the empty space that dimensions vary?

In advance, thanks for your help.

Answer 1

Your function works with both arrays:

In [1]: def foo(V):
   ...:     return V[0]+V[1]
   ...: 
In [2]: foo(np.array([1,3]))
Out[2]: 4
In [3]: foo(np.array([[[1,2],[3,4]], [[5,6],[7,8]]]))
Out[3]: 
array([[ 6,  8],
       [10, 12]])

This answer is just the sum of these two arrays:

In [4]: np.array([[[1,2],[3,4]], [[5,6],[7,8]]])[0]
Out[4]: 
array([[1, 2],
       [3, 4]])
In [5]: np.array([[[1,2],[3,4]], [[5,6],[7,8]]])[1]
Out[5]: 
array([[5, 6],
       [7, 8]])

If you expected something else, you'll have to show us.

As for your second question:

In [6]: t1=np.array([[1,2,3], [4,5,6]])
   ...: t2=np.array([1,2,3])
   ...: t3=np.array([[1,2,3], [4,5,6],5])
   ...: 
In [7]: t1.shape
Out[7]: (2, 3)
In [8]: t2.shape
Out[8]: (3,)
In [9]: t3.shape
Out[9]: (3,)

(3,) is a 1 element tuple. Compare these expressions.

In [11]: (3)
Out[11]: 3
In [12]: (3,)
Out[12]: (3,)

There have been several recent questions about (3,) v (3,1) shape arrays, and np.array([[1,2,3]]) v. np.array([1,2,3]).

t3 is an object dtype array, with 3 elements. The 3 inputs are different length, so it can't create a 2d array. Stay away from this type of array for now. Focus on the simpler arrays.

In [10]: t3
Out[10]: array([[1, 2, 3], [4, 5, 6], 5], dtype=object)
In [13]: t3[0]
Out[13]: [1, 2, 3]
In [14]: t3[2]
Out[14]: 5

Numpy: Why is difference of a (2,1) array and a vertical matrix slice not a (2,1) array

Difference between single and double bracket Numpy array?

=====================

With the nGauss:

In [53]: mu=np.array([0,0])
In [54]: cov=np.eye(2)
In [55]: xx=np.array([[[1,2], [5,6]], [[7,8],[9,0]]])
In [56]: np.apply_along_axis(nGauss, -1, xx, mu, cov)
Out[56]: 
array([[ -1.30642333e-02,  -9.03313360e-15],
       [ -4.61510838e-26,  -4.10103631e-19]])

apply_along_axis iterates on the 1st 2 dim, passing each xx[i,j,:] to nGauss. It's not fast, but is relatively easy to apply.

k = X.shape[0];  # I assume you want
k = X.shape[[1]   # the last dimension
dev = X-mu     # works as long as mu has k terms

this is a scalar:

p1 = np.power( np.power(np.pi * 2, k) , -0.5);

so is

p2 = np.power( np.linalg.det(cov)  , -0.5)

So it comes down to generalizing this expression:

p3 = np.exp( -0.5 * np.dot( np.dot(dev.transpose(), np.linalg.inv(cov)), dev));

In the simple (2,) x case, dev is 1d, and dev.transpose() does nothing.

It's easier to generalize einsum than dot; I think the equivalent is:

p3 = np.einsum('j,j', np.einsum('i,ij', dev, np.linalg.inv(cov)), dev)
p3 = np.exp( -0.5 * p3)

which simplifies to

p3 = np.einsum('i,ij,j', dev, np.linalg.inv(cov), dev)

generalizing to higher dim:

p3 = np.einsum('...i,ij,...j', dev, np.linalg.inv(cov), dev)

So with:

def nGaussA(X, mu, cov):
    # multivariate negative gaussian.    
    # mu is a vector and cov is a covariance matrix.

    k = X.shape[-1];
    dev = X-mu
    p1 = np.power( np.power(np.pi * 2, k) , -0.5);
    p2 = np.power( np.linalg.det(cov)  , -0.5)
    p3 = np.einsum('...i,ij,...j', dev, np.linalg.inv(cov), dev)
    p3 = np.exp( -0.5 * p3)
    return -1.0 * p1 * p2 * p3;

matching earlier values:

In [85]: nGaussA(x,mu,cov)
Out[85]: -0.013064233284684921
In [86]: nGaussA(xx,mu,cov)
Out[86]: 
array([[ -1.30642333e-02,  -9.03313360e-15],
       [ -4.61510838e-26,  -4.10103631e-19]])

So the way to generalize the function is to check each step. If it produces a scalar, keep it. If operates with an x keep it. But if it requires coordinating dimensions with other arrays, use a numpy operation that does that. Often that involves broadcasting. Sometimes it helps to study other numpy functions to see how they generalize (e.g. apply_along_axis, apply_over_axes, cross, etc).

An interactive numpy session is essential; allowing me to try ideas with small sample arrays.

Answer 2

for Q1 you can pack and unpack arguments:

def foo(*args):

    result = []
    for v in args:
        result.append(v[0] + v[1])

    return result

This will allow you pass in as many vector arguments as you want, then iterate over them, returning a list of each result. You can also pack and unpack kwargs with **. More info here:

https://docs.python.org/2/tutorial/controlflow.html#unpacking-argument-lists

Answer 3

For Q1, I'm guessing you want to add the innermost dimensions of your arrays, regardless of how many dimensions the arrays have. The simplest way to do this is to use ellipsis indexing. Here's a detailed example:

>>> a = np.arange(24).reshape((3, 4, 2))
>>> a
array([[[ 0,  1],
        [ 2,  3],
        [ 4,  5],
        [ 6,  7]],

       [[ 8,  9],
        [10, 11],
        [12, 13],
        [14, 15]],

       [[16, 17],
        [18, 19],
        [20, 21],
        [22, 23]]])
>>> a[..., 0]
array([[ 0,  2,  4,  6],
       [ 8, 10, 12, 14],
       [16, 18, 20, 22]])
>>> a[..., 1]
array([[ 1,  3,  5,  7],
       [ 9, 11, 13, 15],
       [17, 19, 21, 23]])
>>> a[..., 0] + a[..., 1]
array([[ 1,  5,  9, 13],
       [17, 21, 25, 29],
       [33, 37, 41, 45]])

This works equally well for a 1D array:

>>> a = np.array([1, 2])
>>> a[..., 0] + a[..., 1]
3

So just define foo as:

def foo(V):
    return V[..., 0] + V[..., 1]

---

For your nGauss function, the simplest solution is to use np.apply_along_axis. For example, you would call it like this:

>>> np.apply_along_axis(nGauss, -1, x1, mu, cov)