headers
Neal Becker | 1 Jan 20:07
Picon

fast iteration (I think I've got it)

This is a c-api question.

I'm trying to get iterators that are both fast and reasonably general.  I
did confirm that iterating using just the general PyArrayIterObject
protocol is not as fast as using c-style pointers for contiguous arrays.

Please confirm if my understanding is correct.  There are 2 cases to
consider.  There are plain old dense arrays, which will be contiguous, and
there are array 'views' or 'slices'.  The views will not be contiguous. 
However, in all cases, the strides between elements in one dimension are
constant.  This last point is key.  As long as the stride in the last
dimension is a constant, a fast iterator is possible.

I put together a small test, which seems to work for the cases I tried.
The idea is to use the general iterator (via PyArray_IterAllButAxis), but
iteration over the last dimension is done with a c-style pointer.

I have tested it with these cases, and the results appear correct:

a = array ((xrange(4),xrange(4,8),xrange(8,12)), int)
b = a[:,::2]
c = a[:,1::2]
d = a[:,-1::-1]
sum3(a)
0 1 2 3 
4 5 6 7 
8 9 10 11 
sum3(b)
0 2 
4 6
(Continue reading)

Permalink | Reply | 3 comments |
headers
Ondrej Certik | 1 Jan 21:07
Picon
Gravatar

Re: planet.scipy.org

On Dec 31, 2007 11:43 PM, Jarrod Millman <millman <at> berkeley.edu> wrote:
> Hey,
>
> I just wanted to announce that we now have a NumPy/SciPy blog
> aggregator thanks to Gaël Varoquaux:  http://planet.scipy.org/

Cool! I was going to write an email that it could be a good idea. Thanks Gael.

> Feel free to contact me if you have a blog that you would like included.

Is it just for anything, or just scipy/numpy related?

I use this blog:

http://ondrejcertik.blogspot.com/

however, since I am new to blogging, I am not sure what the best practise is.

I blog about a lot of things (SymPy, Sage, Debian), but I think just
something is suitable to be on planet.scipy.org.
So I think the best way is to just use tags. So could you please put
this into the planet.conf?

[http://ondrejcertik.blogspot.com/feeds/posts/default/-/scipy]
name = Ond&#345;ej &#268;ert&iacute;k

and I'll be tagging scipy related blogposts with the scipy tag
(currently 0 posts:).

Ondrej
(Continue reading)

Permalink | Reply | 7 comments |
headers
Anne Archibald | 1 Jan 21:34
Picon

Re: fast iteration (I think I've got it)

On 01/01/2008, Neal Becker <ndbecker2 <at> gmail.com> wrote:
> This is a c-api question.
>
> I'm trying to get iterators that are both fast and reasonably general.  I
> did confirm that iterating using just the general PyArrayIterObject
> protocol is not as fast as using c-style pointers for contiguous arrays.

I'd like to point out that "contiguous" can be misleading as it is
used in numpy. An array is flagged contiguous if all the elements are
contiguous *and* they are arranged as a C-ordered multidimensional
array - i.e., the last index changes the fastest as you walk through
the array elements, and the next-to-last chenges the next fastest, and
so on.

> Please confirm if my understanding is correct.  There are 2 cases to
> consider.  There are plain old dense arrays, which will be contiguous, and
> there are array 'views' or 'slices'.  The views will not be contiguous.
> However, in all cases, the strides between elements in one dimension are
> constant.  This last point is key.  As long as the stride in the last
> dimension is a constant, a fast iterator is possible.
>
> I put together a small test, which seems to work for the cases I tried.
> The idea is to use the general iterator (via PyArray_IterAllButAxis), but
> iteration over the last dimension is done with a c-style pointer.

This is not an unreasonable approach, but it can suffer badly if the
last dimension is "small", for example if you have an array of shape
(1000,1000,2). Such arrays can arise for example from working with
complex numbers as pairs of real numbers.
(Continue reading)

Permalink | Reply | 1 comment |
headers
Neal Becker | 1 Jan 22:22
Picon

Re: fast iteration (I think I've got it)

Thank you for the response.  I'm afraid I haven't explained what I'm doing.

I have a lot of c++ code written in a generic c++ interface style.  (The
code is for signal processing, but that is irrelevant). 

The code is generic for data types as well as container types.

To accomplish this, the interface uses the boost::range interface concept. 
An example of using this interface:

template<typename in_t>
void F (in_t const& in) {
  typename boost::range_const_iterator<in_t>::type i = boost::begin (in);
  for (; i != boost::end (in); ++i)
    do_something_with (*i);
}

In the above, in_t is a container type.  It could be std::vector<int>, for
example.

The concept has:
  range_iterator<in_t>::type and range_const_iterator<in_t>::type are the
iterator types for the range
  begin(in) gives and iterator pointing to the beginning
  ++i increments the iterator
  *i derefs the iterator

This allows writing functions that work with different container types.  For
example, std::vector, boost::ublas::vector.
(Continue reading)

Permalink | Reply | 0 comments |
headers
Jarrod Millman | 1 Jan 22:31
Picon
Favicon
Gravatar

Re: planet.scipy.org

On Jan 1, 2008 12:07 PM, Ondrej Certik <ondrej <at> certik.cz> wrote:
> So could you please put this into the planet.conf?
>
> [http://ondrejcertik.blogspot.com/feeds/posts/default/-/scipy]
> name = Ond&#345;ej &#268;ert&iacute;k

Done.
Thanks,

--

-- 
Jarrod Millman
Computational Infrastructure for Research Labs
10 Giannini Hall, UC Berkeley
phone: 510.643.4014
http://cirl.berkeley.edu/
Permalink | Reply | 6 comments |
headers
Travis E. Oliphant | 1 Jan 22:37

Re: fast iteration (I think I've got it)

Neal Becker wrote:
> This is a c-api question.
>
> I'm trying to get iterators that are both fast and reasonably general.  I
> did confirm that iterating using just the general PyArrayIterObject
> protocol is not as fast as using c-style pointers for contiguous arrays.
>
> Please confirm if my understanding is correct.  There are 2 cases to
> consider.  There are plain old dense arrays, which will be contiguous, and
> there are array 'views' or 'slices'.  The views will not be contiguous. 
> However, in all cases, the strides between elements in one dimension are
> constant.  This last point is key.  As long as the stride in the last
> dimension is a constant, a fast iterator is possible.
>   

Your understanding is correct.  You can use the iterator over all but 
one dimension.  This is a common paradigm in NumPy (use an iterator over 
all but one dimension (usually chosen to be the one with the smallest 
striding) and then use strided access over the final dimension (see all 
the ufuncs for examples).

IterAllButAxis takes a pointer to an integer.  If that integer is 
negative on input, then IterAllButAxis will choose the axis for you to 
be the one with the smallest stride,  then it will set the variable to 
the axis that was chosen. 

This is just what ufuncs do and is a pretty good idea generally to 
minimize the over-head of the final inner loop.

-Travis O.
(Continue reading)

Permalink | Reply | 0 comments |
headers
Jarrod Millman | 2 Jan 02:14
Picon
Favicon
Gravatar

Re: [OT] Which version of svn does svn.scipy.org run ?

On Dec 17, 2007 11:23 PM, David Cournapeau <david <at> ar.media.kyoto-u.ac.jp> wrote:
>     Is it possible to know which version of subversion is used by the
> svn.scipy.org server ? I am trying to use svnsync on it, without any
> success, and wanted to know if this was coming from svn.scipy.org or
> from something else ?

It is running subversion 1.2.3-2.1.  In the near future, the server
will be replaced with newer hardware and updated software.

--

-- 
Jarrod Millman
Computational Infrastructure for Research Labs
10 Giannini Hall, UC Berkeley
phone: 510.643.4014
http://cirl.berkeley.edu/
Permalink | Reply | 1 comment |
headers
Barry Wark | 2 Jan 03:18
Picon
Gravatar

Re: planet.scipy.org

Jarrod,

I'd like to submit my blog for inclusion too...

http://softwareforscientists.blogspot.com/search/label/scipy

Again, not too many posts yet, but things will be ramping up.

On Jan 1, 2008 1:31 PM, Jarrod Millman <millman <at> berkeley.edu> wrote:
> On Jan 1, 2008 12:07 PM, Ondrej Certik <ondrej <at> certik.cz> wrote:
> > So could you please put this into the planet.conf?
> >
> > [http://ondrejcertik.blogspot.com/feeds/posts/default/-/scipy]
> > name = Ond&#345;ej &#268;ert&iacute;k
>
> Done.
> Thanks,
>
> --
> Jarrod Millman
> Computational Infrastructure for Research Labs
> 10 Giannini Hall, UC Berkeley
> phone: 510.643.4014
> http://cirl.berkeley.edu/
> _______________________________________________
>
> Numpy-discussion mailing list
> Numpy-discussion <at> scipy.org
> http://projects.scipy.org/mailman/listinfo/numpy-discussion
>
(Continue reading)

Permalink | Reply | 5 comments |
headers
Jarrod Millman | 2 Jan 03:35
Picon
Favicon
Gravatar

Re: planet.scipy.org

On Jan 1, 2008 6:18 PM, Barry Wark <barrywark <at> gmail.com> wrote:
> I'd like to submit my blog for inclusion too...
>
> http://softwareforscientists.blogspot.com/search/label/scipy

Hmm.  I added you, but while I was adding you it looks like your blog
disappeared.
Thanks,

--

-- 
Jarrod Millman
Computational Infrastructure for Research Labs
10 Giannini Hall, UC Berkeley
phone: 510.643.4014
http://cirl.berkeley.edu/
Permalink | Reply | 4 comments |
headers
Neal Becker | 2 Jan 03:49
Picon

need inverse of PyArray_ITER_GOTO1D

I need a function that is the inverse of PyArray_ITER_GOTO1D.  A bit of
guesswork brought me to this, can anyone please tell me if it looks
correct?

inline npy_intp as_linear (PyArrayIterObject const* it) {
  if (it->nd_m1 == 0)
    return (it->dataptr - it->ao->data)/it->strides[0];
  else if (it->contiguous)
    return (it->dataptr - it->ao->data)/it->ao->descr->elsize;
  else {
    npy_intp loc = 0;
    npy_intp offset = (it->dataptr - it->ao->data)/it->ao->descr->elsize;
    for (int i = 0; i <= it->nd_m1; ++i) {
      loc += offset / it->factors[i];
      offset %= it->factors[i];
    }
    return loc;
  }
}
Permalink | Reply | 1 comment |

Gmane