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:
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 |
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:
_______________________________________________ SciPy-User mailing list [hidden email] http://mail.scipy.org/mailman/listinfo/scipy-user |
Dan, Thanks for your answer. 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:
_______________________________________________ SciPy-User mailing list [hidden email] http://mail.scipy.org/mailman/listinfo/scipy-user |
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 |
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 |
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 |
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:
_______________________________________________ SciPy-User mailing list [hidden email] http://mail.scipy.org/mailman/listinfo/scipy-user |
Free forum by Nabble | Edit this page |