Numpy/Scipy: Avoiding nested loops to operate on matrix-valued images

classic Classic list List threaded Threaded
7 messages Options
Reply | Threaded
Open this post in threaded view
|

Numpy/Scipy: Avoiding nested loops to operate on matrix-valued images

tyldurd

Hello,

I am a beginner at python and numpy and I need to compute the matrix logarithm for each "pixel" (i.e. x,y position) of a matrix-valued image of dimension MxNx3x3. 3x3 is the dimensions of the matrix at each pixel.

The function I have written so far is the following:

def logm_img(im):
   
from scipy import linalg
    dimx
= im.shape[0]
    dimy
= im.shape[1]
    res
= zeros_like(im)
   
for x in range(dimx):
       
for y in range(dimy):
            res
[x, y, :, :] = linalg.logm(asmatrix(im[x,y,:,:]))
   
return res

Is it ok? Is there a way to avoid the two nested loops ?


_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: Numpy/Scipy: Avoiding nested loops to operate on matrix-valued images

Dan Lussier
Have you tried numpy.frompyfunc?


With this approach you may be able create a function which acts elementwise over your array to compute the matrix logarithm at each entry using Numpy's ufuncs.  This would avoid the explicit iteration over the array using the for loops.

As a rough outline try:

from scipy import linalg
import numpy as np

# Assume im is the container array containing a 3x3 matrix at each pixel.

# Composite function so get matrix log of array A as a matrix in one step
def log_matrix(A):
    return linalg.logm(np.asmatrix(A))
   

# Creating function to operate over container array.  Takes one argument and returns the result.
log_ufunc = np.frompyfunc(log_matrix, 1, 1)

# Using log_ufunc on container array, im
res = log_ufunc(im)

Dan
   

On 2012-03-15, at 1:59 AM, tyldurd wrote:

Hello,

I am a beginner at python and numpy and I need to compute the matrix logarithm for each "pixel" (i.e. x,y position) of a matrix-valued image of dimension MxNx3x3. 3x3 is the dimensions of the matrix at each pixel.

The function I have written so far is the following:

def logm_img(im):
   
from scipy import linalg
    dimx
= im.shape[0]
    dimy
= im.shape[1]
    res
= zeros_like(im)
   
for x in range(dimx):
       
for y in range(dimy):
            res
[x, y, :, :] = linalg.logm(asmatrix(im[x,y,:,:]))
   
return res

Is it ok? Is there a way to avoid the two nested loops ?

_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user


_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: Numpy/Scipy: Avoiding nested loops to operate on matrix-valued images

tyldurd
Dan,

Thanks for your answer.

However, this solution does not work for me. 

First, it returns an array with dtype=object which is not the original type of the data.  Besides, the values in the array are not equal to the ones given by the 'traditional' nested loops. I think the problem comes from the fact that ufuncs are functions that act over each element of an array, not over slices.

I have done a lot of research on this topic but it seems it is not feasible in terms of slicing or vectorizing. The only solution I found would be with generalized ufuncs but from what I understand, they require to write C code, which I would like to avoid :-)

Therefore, I am going to stick to nested loops at least for now.

Regards,

Olivier

On Thursday, March 15, 2012 9:23:49 PM UTC+1, Dan Lussier wrote:
Have you tried numpy.frompyfunc?


With this approach you may be able create a function which acts elementwise over your array to compute the matrix logarithm at each entry using Numpy's ufuncs.  This would avoid the explicit iteration over the array using the for loops.

As a rough outline try:

from scipy import linalg
import numpy as np

# Assume im is the container array containing a 3x3 matrix at each pixel.

# Composite function so get matrix log of array A as a matrix in one step
def log_matrix(A):
    return linalg.logm(np.asmatrix(A))
   

# Creating function to operate over container array.  Takes one argument and returns the result.
log_ufunc = np.frompyfunc(log_matrix, 1, 1)

# Using log_ufunc on container array, im
res = log_ufunc(im)

Dan
   

On 2012-03-15, at 1:59 AM, tyldurd wrote:

Hello,

I am a beginner at python and numpy and I need to compute the matrix logarithm for each "pixel" (i.e. x,y position) of a matrix-valued image of dimension MxNx3x3. 3x3 is the dimensions of the matrix at each pixel.

The function I have written so far is the following:

def logm_img(im):
   
from scipy import linalg
    dimx
= im.shape[0]
    dimy
= im.shape[1]
    res
= zeros_like(im)
   
for x in range(dimx):
       
for y in range(dimy):
            res
[x, y, :, :] = linalg.logm(asmatrix(im[x,y,:,:]))
   
return res

Is it ok? Is there a way to avoid the two nested loops ?

_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user


_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: Numpy/Scipy: Avoiding nested loops to operate on matrix-valued images

Nathaniel Smith
On Mon, Mar 19, 2012 at 9:23 AM, tyldurd <[hidden email]> wrote:
> I have done a lot of research on this topic but it seems it is not feasible
> in terms of slicing or vectorizing. The only solution I found would be with
> generalized ufuncs but from what I understand, they require to write C code,
> which I would like to avoid :-)

I think the idea of generalized ufuncs is that linalg.logm should be
written as a generalized ufunc already out of the box, and then this
would be straightforward. However: (1) it isn't, and (2) even if it
were, I'm having trouble understanding from the available docs how you
would actually use it -- maybe calling logm would just work for your
case, but there don't seem to be any examples available of how it
chooses which dimensions to apply to. (Are there any generalized
ufuncs actually defined in the standard packages? For instance, is
np.dot implemented as a generalized ufunc? Should it be?)

> Therefore, I am going to stick to nested loops at least for now.

That seems like the best option to me. Nothing immoral about using a
loop when that's what you need :-).

-- Nathaniel
_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: Numpy/Scipy: Avoiding nested loops to operate on matrix-valued images

josef.pktd
On Tue, Mar 20, 2012 at 11:33 AM, Nathaniel Smith <[hidden email]> wrote:

> On Mon, Mar 19, 2012 at 9:23 AM, tyldurd <[hidden email]> wrote:
>> I have done a lot of research on this topic but it seems it is not feasible
>> in terms of slicing or vectorizing. The only solution I found would be with
>> generalized ufuncs but from what I understand, they require to write C code,
>> which I would like to avoid :-)
>
> I think the idea of generalized ufuncs is that linalg.logm should be
> written as a generalized ufunc already out of the box, and then this
> would be straightforward. However: (1) it isn't, and (2) even if it
> were, I'm having trouble understanding from the available docs how you
> would actually use it -- maybe calling logm would just work for your
> case, but there don't seem to be any examples available of how it
> chooses which dimensions to apply to. (Are there any generalized
> ufuncs actually defined in the standard packages? For instance, is
> np.dot implemented as a generalized ufunc? Should it be?)

only in a test case, AFAIK
from numpy.core.umath_tests import matrix_multiply

Josef

>
>> Therefore, I am going to stick to nested loops at least for now.
>
> That seems like the best option to me. Nothing immoral about using a
> loop when that's what you need :-).
>
> -- Nathaniel
> _______________________________________________
> SciPy-User mailing list
> [hidden email]
> http://mail.scipy.org/mailman/listinfo/scipy-user
_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: Numpy/Scipy: Avoiding nested loops to operate on matrix-valued images

David Warde-Farley-3
In reply to this post by Nathaniel Smith
On 2012-03-20, at 11:33 AM, Nathaniel Smith wrote:

>  (Are there any generalized
> ufuncs actually defined in the standard packages? For instance, is
> np.dot implemented as a generalized ufunc? Should it be?)

Ideally, so long as it still made use of BLAS for the actual matrix products.

I tried my hand at implementing a gufunc for log(sum(exp(...))), with the sum being the "generalized" part. Did not have much luck...

David
_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: Numpy/Scipy: Avoiding nested loops to operate on matrix-valued images

mdekauwe
In reply to this post by tyldurd
I didn't quite follow exactly what you were doing, but someone previously showed me how to avoid inner loops and so perhaps this will help?

Instead of...

tmp = np.arange(500000).reshape(1000,500)
nrows, ncols = tmp.shape[0], tmp.shape[1]
out = np.zeros((nrows, ncols))
for i in xrange(nrows):
    for j in xrange(ncols):
        out[i,j] = tmp[i,j] * 3.0

you might try...

tmp = np.arange(500000).reshape(1000,500)
nrows, ncols = tmp.shape[0], tmp.shape[1]
out = np.zeros((nrows, ncols))
r = np.arange(nrows)
c = np.arange(ncols)
out[r[:,None],c] = tmp[r[:,None],c] * 3.0

Assuming your arrays are large you would get a speed bump


On Thursday, March 15, 2012 7:59:28 PM UTC+11, tyldurd wrote:

Hello,

I am a beginner at python and numpy and I need to compute the matrix logarithm for each "pixel" (i.e. x,y position) of a matrix-valued image of dimension MxNx3x3. 3x3 is the dimensions of the matrix at each pixel.

The function I have written so far is the following:

def logm_img(im):
   
from scipy import linalg
    dimx
= im.shape[0]
    dimy
= im.shape[1]
    res
= zeros_like(im)
   
for x in range(dimx):
       
for y in range(dimy):
            res
[x, y, :, :] = linalg.logm(asmatrix(im[x,y,:,:]))
   
return res

Is it ok? Is there a way to avoid the two nested loops ?


_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user