Making A Function That Can Take Arguments In Various Shapes

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

Solution 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]: 3In[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]: t3Out[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 wantk = 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:

defnGaussA(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.

Solution 2:

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:

deffoo(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)

Solution 3:

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