This is probably an astonishingly simple question, but I've been struggling for some time with it. I am trying to work with weave for the first time, and it is becoming clear that I don't understand how 2D arrays get passed. The example below illustrates the problem I'm having... the 2D indexing doesn't seem to be working right at all. BTW the link to the weave documentation is a dead link. Is there any documentation beyond a couple of short examples? Every promising google link turns up dead. #!/usr/bin/env python from numpy import array, zeros, abs, ones, arange, ravel, nan, empty, isnan, histogram, hstack, multiply, log, shape, triu import scipy.weave as weave from scipy.weave import converters def Matrix(width, minval, maxval, exponent, type="HorseSaddle"): halfSqrt2 = 0.707106781185 sin45 = halfSqrt2 cos45 = halfSqrt2 Matrix = zeros((width, width), dtype=float) for j in range(0,width): jScaled = j / (halfSqrt2 * (width-1)) jTranslated = jScaled - halfSqrt2 jsin45 = jTranslated * sin45 jcos45 = jsin45 for i in range(j,width): iScaled = i / (halfSqrt2 * (width-1)); iTranslated = iScaled - halfSqrt2; iRotated = iTranslated * cos45 + jsin45; jRotated = - iTranslated * sin45 + jcos45; iRotated = abs(iRotated); jRotated = abs(jRotated); Matrix[i, j] = round(maxval * iRotated**exponent + \ minval * jRotated**exponent, 2) Matrix[j,i] = Matrix[i,j] return Matrix ################################################################################ # algorithm - WEAVE version ################################################################################ def Fast(s1, s2, wt1, wt2, M): seq1 = s1[0] seq2 = s2[0] t1 = s1[1] t2 = s2[1] code = ''' for (unsigned int idx2 = 0; idx2 < 10; ++idx2) { printf("\\n"); for (unsigned int idx1 = 0; idx1 < 10; ++idx1) { printf("%6f ", M[idx1, idx2]); }// for whole s1 }// for whole s2 printf("\\n"); ''' weave.inline(code, ["seq1", "seq2", "t1", "t2", "M" ]) ################################################################################ # test section ################################################################################ if __name__ == '__main__' : attr1 = array([1,2,2,3,2,4,3,5,4,3,3, 3,4,6,1,0,2,4,2,5,2,3]) attr2 = array([1,2,2,3,2,5,3,5,4,3,2,6, 3,4,6,1,0,1,4,2,5,2,3]) t = arange(1, len(attr1)+1, dtype=float) wt1 = ones(len(attr1), dtype=float) t2 = arange(1, len(attr2)+1, dtype=float) wt2 = ones(len(attr2), dtype=float) seq1 = [attr1, t] seq2 = [attr2, t2] M = Matrix(7,-10,10,2.) print M s1s2, s2s1 = Fast(seq1, seq2, wt1, wt2, M ) -- ----------------------------------------------------------------------- | Alan K. Jackson | To see a World in a Grain of Sand | | alan@ajackson.org | And a Heaven in a Wild Flower, | | www.ajackson.org | Hold Infinity in the palm of your hand | | Houston, Texas | And Eternity in an hour. - Blake | -----------------------------------------------------------------------
On Sat, Nov 17, 2007 at 03:30:10PM -0600, Alan Jackson wrote:
This is probably an astonishingly simple question, but I've been struggling for some time with it.
I am trying to work with weave for the first time, and it is becoming clear that I don't understand how 2D arrays get passed.
Me neither. This is why I use the blitz converter. I modified your example to use the vlitz converters, it is just so much nicer: ++++++++++++++++++++++++++++++++++++++++++++ def Fast(s1, s2, wt1, wt2, M): seq1 = s1[0] seq2 = s2[0] t1 = s1[1] t2 = s2[1] code = ''' for (int idx2 = 0; idx2 < 10; ++idx2) { printf("\\n"); for (int idx1 = 0; idx1 < 10; ++idx1) { printf("%6f ", M(idx1, idx2)); }// for whole s1 }// for whole s2 printf("\\n"); ''' weave.inline(code, ["seq1", "seq2", "t1", "t2", "M" ], type_converters=converters.blitz, ) ++++++++++++++++++++++++++++++++++++++++++++
BTW the link to the weave documentation is a dead link. Is there any documentation beyond a couple of short examples? Every promising google link turns up dead.
There have been a few threads on this mailing list. I also think http://scipy.org/PerformancePython#head-a3f4dd816378d3ba4cbdd3d23dc98529e8ad... can be useful. Gaël
Gael, have you tried weave with openmp? Cheers, William On Nov 17, 2007 4:46 PM, Gael Varoquaux <gael.varoquaux@normalesup.org> wrote:
On Sat, Nov 17, 2007 at 03:30:10PM -0600, Alan Jackson wrote:
This is probably an astonishingly simple question, but I've been struggling for some time with it.
I am trying to work with weave for the first time, and it is becoming clear that I don't understand how 2D arrays get passed.
Me neither. This is why I use the blitz converter. I modified your example to use the vlitz converters, it is just so much nicer:
++++++++++++++++++++++++++++++++++++++++++++ def Fast(s1, s2, wt1, wt2, M):
seq1 = s1[0] seq2 = s2[0] t1 = s1[1] t2 = s2[1]
code = '''
for (int idx2 = 0; idx2 < 10; ++idx2) { printf("\\n"); for (int idx1 = 0; idx1 < 10; ++idx1) { printf("%6f ", M(idx1, idx2));
}// for whole s1 }// for whole s2 printf("\\n"); '''
weave.inline(code, ["seq1", "seq2", "t1", "t2", "M" ], type_converters=converters.blitz, ) ++++++++++++++++++++++++++++++++++++++++++++
BTW the link to the weave documentation is a dead link. Is there any documentation beyond a couple of short examples? Every promising google link turns up dead.
There have been a few threads on this mailing list. I also think http://scipy.org/PerformancePython#head-a3f4dd816378d3ba4cbdd3d23dc98529e8ad... can be useful.
Gaël
_______________________________________________ SciPy-user mailing list SciPy-user@scipy.org http://projects.scipy.org/mailman/listinfo/scipy-user
On Sat, Nov 17, 2007 at 05:16:55PM -0500, william ratcliff wrote:
Gael, have you tried weave with openmp?
No, and I would be interested to find out if it works or not. But if it doesn't, I don't see why ctypes would not (maybe I am missing some thing, though, my knowledge of dynamical loading is very thin). Gaël
On Sat, 17 Nov 2007 22:46:59 +0100 Gael Varoquaux <gael.varoquaux@normalesup.org> wrote:
On Sat, Nov 17, 2007 at 03:30:10PM -0600, Alan Jackson wrote:
This is probably an astonishingly simple question, but I've been struggling for some time with it.
I am trying to work with weave for the first time, and it is becoming clear that I don't understand how 2D arrays get passed.
Me neither. This is why I use the blitz converter. I modified your example to use the vlitz converters, it is just so much nicer:
Well, with my *real* code, blitz gave me several pages of errors - I tried treating the 2D arrays as a 1D array, doing all the index arithmetic myself, and that seems to work. Which strikes me as odd... -- ----------------------------------------------------------------------- | Alan K. Jackson | To see a World in a Grain of Sand | | alan@ajackson.org | And a Heaven in a Wild Flower, | | www.ajackson.org | Hold Infinity in the palm of your hand | | Houston, Texas | And Eternity in an hour. - Blake | -----------------------------------------------------------------------
On Sat, Nov 17, 2007 at 04:44:52PM -0600, Alan Jackson wrote:
On Sat, 17 Nov 2007 22:46:59 +0100 Well, with my *real* code, blitz gave me several pages of errors - I tried treating the 2D arrays as a 1D array, doing all the index arithmetic myself, and that seems to work. Which strikes me as odd...
With Blitz, you are allowed M(i,j) indexing. Otherwise, only the memory address of the data is passed and, as you noticed, you must calculate the offset yourself, i.e. M[i*nr_cols + j]. I modified your code to work with ctypes, which you may find interesting. Use python setup.py build_ext -i to build the C file. Then run alan.py. Cheers Stéfan
Stefan van der Walt wrote:
On Sat, Nov 17, 2007 at 04:44:52PM -0600, Alan Jackson wrote:
On Sat, 17 Nov 2007 22:46:59 +0100 Well, with my *real* code, blitz gave me several pages of errors - I tried treating the 2D arrays as a 1D array, doing all the index arithmetic myself, and that seems to work. Which strikes me as odd...
With Blitz, you are allowed M(i,j) indexing. Otherwise, only the memory address of the data is passed and, as you noticed, you must calculate the offset yourself, i.e. M[i*nr_cols + j].
Travis added a few macros to automate this calculation (for non-blitz cases). If you passed in 'a' then A1(i), A2(i, j) A3(i, j, k) and A4(i, j, k, l) are all defined. Here are the macros: #define A1(i) (*((long*)(a_array->data + (i)*Sa[0]))) #define A2(i,j) (*((long*)(a_array->data + (i)*Sa[0] + (j)*Sa[1]))) #define A3(i,j,k) (*((long*)(a_array->data + (i)*Sa[0] + (j)*Sa[1] + (k)*Sa[2]))) #define A4(i,j,k,l) (*((long*)(a_array->data + (i)*Sa[0] + (j)*Sa[1] + (k)*Sa[2] + (l)*Sa[3]))) cheers, prabhu
I had good luck with it using blitz type converters (though if there's a bug in the C code, it's painful to find ;>). Did you look at the array3d example that came with the code? So far, I've been able to pass 1,2, and 3d arrays successfully. But, I've been cheating and creating the arrays that I want returned in numpy before calling the routine in weave and passing them in as arguments. I modify them within weave, and then check the results when weave terminates to make sure they have the values I think they should have in numpy. Then I didn't have to work on figuring out return types, allocating memory,etc. I also use blitz type converters. What kind of errors are you getting? When I was writing my first program using weave, I got tons of errors--even if it was something as minor as forgetting a ";" in one of the lines of code--one helpful thing to do is to find the temporary directory where weave stores the c++ versions of your code and deleting it--there should be two directories. However, if it's not a simple question of stale code lying around, then wade through the pages of output and try to find out whether there is a simple bug in your C code. The next possible source of error could be a compiler issue--for example, I normally use an older version of gcc (I forget, something like 3.4) and upgraded to version 4.x (one of the most recent) and found that it no longer worked. Is the real code involved something simple enough that you can post it? Cheers, William On Nov 17, 2007 5:44 PM, Alan Jackson <alan@ajackson.org> wrote:
On Sat, 17 Nov 2007 22:46:59 +0100 Gael Varoquaux <gael.varoquaux@normalesup.org> wrote:
On Sat, Nov 17, 2007 at 03:30:10PM -0600, Alan Jackson wrote:
This is probably an astonishingly simple question, but I've been struggling for some time with it.
I am trying to work with weave for the first time, and it is becoming clear that I don't understand how 2D arrays get passed.
Me neither. This is why I use the blitz converter. I modified your example to use the vlitz converters, it is just so much nicer:
Well, with my *real* code, blitz gave me several pages of errors - I tried treating the 2D arrays as a 1D array, doing all the index arithmetic myself, and that seems to work. Which strikes me as odd...
-- ----------------------------------------------------------------------- | Alan K. Jackson | To see a World in a Grain of Sand | | alan@ajackson.org | And a Heaven in a Wild Flower, | | www.ajackson.org | Hold Infinity in the palm of your hand | | Houston, Texas | And Eternity in an hour. - Blake | ----------------------------------------------------------------------- _______________________________________________
SciPy-user mailing list SciPy-user@scipy.org http://projects.scipy.org/mailman/listinfo/scipy-user
Thanks all for your help. The path I chose was to redo my indexing - that was the path of least resistance. I'm not sure where the blitz errors came from - the same code compiled clean with inline, so that is a mystery (I did modify the arrays for blitz before compiling). I also created all the arrays I wanted returned in numpy first - which makes for a really long function call, but it works and it is not too awkward. I got a speed up of x37 on my real test. On Sat, 17 Nov 2007 18:16:26 -0500 "william ratcliff" <william.ratcliff@gmail.com> wrote:
I had good luck with it using blitz type converters (though if there's a bug in the C code, it's painful to find ;>). Did you look at the array3d example that came with the code? So far, I've been able to pass 1,2, and 3d arrays successfully. But, I've been cheating and creating the arrays that I want returned in numpy before calling the routine in weave and passing them in as arguments. I modify them within weave, and then check the results when weave terminates to make sure they have the values I think they should have in numpy. Then I didn't have to work on figuring out return types, allocating memory,etc. I also use blitz type converters.
What kind of errors are you getting? When I was writing my first program using weave, I got tons of errors--even if it was something as minor as forgetting a ";" in one of the lines of code--one helpful thing to do is to find the temporary directory where weave stores the c++ versions of your code and deleting it--there should be two directories. However, if it's not a simple question of stale code lying around, then wade through the pages of output and try to find out whether there is a simple bug in your C code. The next possible source of error could be a compiler issue--for example, I normally use an older version of gcc (I forget, something like 3.4) and upgraded to version 4.x (one of the most recent) and found that it no longer worked. Is the real code involved something simple enough that you can post it?
Cheers, William
On Nov 17, 2007 5:44 PM, Alan Jackson <alan@ajackson.org> wrote:
On Sat, 17 Nov 2007 22:46:59 +0100 Gael Varoquaux <gael.varoquaux@normalesup.org> wrote:
On Sat, Nov 17, 2007 at 03:30:10PM -0600, Alan Jackson wrote:
This is probably an astonishingly simple question, but I've been struggling for some time with it.
I am trying to work with weave for the first time, and it is becoming clear that I don't understand how 2D arrays get passed.
Me neither. This is why I use the blitz converter. I modified your example to use the vlitz converters, it is just so much nicer:
Well, with my *real* code, blitz gave me several pages of errors - I tried treating the 2D arrays as a 1D array, doing all the index arithmetic myself, and that seems to work. Which strikes me as odd...
-- ----------------------------------------------------------------------- | Alan K. Jackson | To see a World in a Grain of Sand | | alan@ajackson.org | And a Heaven in a Wild Flower, | | www.ajackson.org | Hold Infinity in the palm of your hand | | Houston, Texas | And Eternity in an hour. - Blake | ----------------------------------------------------------------------- _______________________________________________
SciPy-user mailing list SciPy-user@scipy.org http://projects.scipy.org/mailman/listinfo/scipy-user
_______________________________________________ SciPy-user mailing list SciPy-user@scipy.org http://projects.scipy.org/mailman/listinfo/scipy-user
-- ----------------------------------------------------------------------- | Alan K. Jackson | To see a World in a Grain of Sand | | alan@ajackson.org | And a Heaven in a Wild Flower, | | www.ajackson.org | Hold Infinity in the palm of your hand | | Houston, Texas | And Eternity in an hour. - Blake | -----------------------------------------------------------------------
On Sat, Nov 17, 2007 at 05:30:18PM -0600, Alan Jackson wrote:
I'm not sure where the blitz errors came from - the same code compiled clean with inline, so that is a mystery (I did modify the arrays for blitz before compiling).
Just checking: you did replace the "M[x, y]" by "M(x, y)" in all your code?
I also created all the arrays I wanted returned in numpy first - which makes for a really long function call, but it works and it is not too awkward.
I always do that. It removes all memory management from the C level. Tht's one less hurdle. I usual wrap my function in a another one, to hide this ugly function call, though.
I got a speed up of x37 on my real test.
Nice ! Gaël
On Sun, 18 Nov 2007 00:36:06 +0100 Gael Varoquaux <gael.varoquaux@normalesup.org> wrote:
On Sat, Nov 17, 2007 at 05:30:18PM -0600, Alan Jackson wrote:
I'm not sure where the blitz errors came from - the same code compiled clean with inline, so that is a mystery (I did modify the arrays for blitz before compiling).
Just checking: you did replace the "M[x, y]" by "M(x, y)" in all your code?
Yep...
I also created all the arrays I wanted returned in numpy first - which makes for a really long function call, but it works and it is not too awkward.
I always do that. It removes all memory management from the C level. Tht's one less hurdle. I usual wrap my function in a another one, to hide this ugly function call, though.
I got a speed up of x37 on my real test.
Nice !
Gaël _______________________________________________ SciPy-user mailing list SciPy-user@scipy.org http://projects.scipy.org/mailman/listinfo/scipy-user
-- ----------------------------------------------------------------------- | Alan K. Jackson | To see a World in a Grain of Sand | | alan@ajackson.org | And a Heaven in a Wild Flower, | | www.ajackson.org | Hold Infinity in the palm of your hand | | Houston, Texas | And Eternity in an hour. - Blake | -----------------------------------------------------------------------
On Sat, Nov 17, 2007 at 04:44:52PM -0600, Alan Jackson wrote:
Well, with my *real* code, blitz gave me several pages of errors
If you are doing complicated things, you should really use another technique to call C from Python. Ctypes is pretty easy to use (http://www.scipy.org/Cookbook/Ctypes). Debug is also much easier, because it is easier to understand the compiler's error. Gaël
I'll have to look into it... On Sun, 18 Nov 2007 00:33:06 +0100 Gael Varoquaux <gael.varoquaux@normalesup.org> wrote:
On Sat, Nov 17, 2007 at 04:44:52PM -0600, Alan Jackson wrote:
Well, with my *real* code, blitz gave me several pages of errors
If you are doing complicated things, you should really use another technique to call C from Python. Ctypes is pretty easy to use (http://www.scipy.org/Cookbook/Ctypes). Debug is also much easier, because it is easier to understand the compiler's error.
Gaël _______________________________________________ SciPy-user mailing list SciPy-user@scipy.org http://projects.scipy.org/mailman/listinfo/scipy-user
-- ----------------------------------------------------------------------- | Alan K. Jackson | To see a World in a Grain of Sand | | alan@ajackson.org | And a Heaven in a Wild Flower, | | www.ajackson.org | Hold Infinity in the palm of your hand | | Houston, Texas | And Eternity in an hour. - Blake | -----------------------------------------------------------------------
participants (5)
-
Alan Jackson -
Gael Varoquaux -
Prabhu Ramachandran -
Stefan van der Walt -
william ratcliff