Say I have a matrix m whose diagonal values I want to replace with the contents of vector new_diagonal. Right now, I'm using m -= diag(m) m += diag(new_diagonal) but that seems pretty wasteful (especially if diag() creates new arrays, which I suspect it does). Is there a standard (or any efficient) way to assign the diagonal values directly without an explicit for loop? Googling for "site:www.scipy.org assign diagonal" only returned results for sparse matrices, which is not the case here - these are plain numpy arrays. Thanks in advance, -carlos
Hello, m -= diag(m) subtracts diag(m) from each row of m for me, if m is a numpy array, so your method could be far worse than wasteful. There may be an easier way to implement, but at its root it will be similar to the simple for loop: for i in range(len(new_diagonal)): m[i][i] = new_diagonal[i] assuming that new_diagonal is already a 1D array. If not, try: new_diagonal = diag(new_diagonal) first then run. If your dimensions are (very) large, you may have increased performance with for i in xrange(len(new_diagonal)): Hope this is helpful, Mclean On Wed, 30 Jan 2008, Carlos Scheidegger wrote:
Say I have a matrix m whose diagonal values I want to replace with the contents of vector new_diagonal. Right now, I'm using
m -= diag(m) m += diag(new_diagonal)
but that seems pretty wasteful (especially if diag() creates new arrays, which I suspect it does). Is there a standard (or any efficient) way to assign the diagonal values directly without an explicit for loop? Googling for "site:www.scipy.org assign diagonal" only returned results for sparse matrices, which is not the case here - these are plain numpy arrays.
Thanks in advance, -carlos _______________________________________________ SciPy-user mailing list SciPy-user@scipy.org http://projects.scipy.org/mailman/listinfo/scipy-user
Try the function below (from pyGAUSS). hth, Alan Isaac #diagrv: insert v as diagonal of matrix x (2D only!) def diagrv(x,v,copy=True): assert(len(x.shape)==2), "For 2-d arrays only." x = numpy.matrix( x, copy=copy ) stride = 1 + x.shape[1] x.flat[ slice(0,None,stride) ] = v return x
On 30-Jan-08, at 10:08 PM, Alan G Isaac wrote:
#diagrv: insert v as diagonal of matrix x (2D only!) def diagrv(x,v,copy=True): assert(len(x.shape)==2), "For 2-d arrays only." x = numpy.matrix( x, copy=copy ) stride = 1 + x.shape[1] x.flat[ slice(0,None,stride) ] = v return x
Ooh that is clever. This should really go on the Cookbook page. Perhaps more generally, for in place modification: def setdiag(m, d): assert(len(x.shape) == 2) stride = 1 + x.shape[1] m.flat[slice(0,None,stride)] = d Regards, David
On Jan 30, 2008 7:44 PM, David Warde-Farley <dwf@cs.toronto.edu> wrote:
On 30-Jan-08, at 10:08 PM, Alan G Isaac wrote:
#diagrv: insert v as diagonal of matrix x (2D only!) def diagrv(x,v,copy=True): assert(len(x.shape)==2), "For 2-d arrays only." x = numpy.matrix( x, copy=copy ) stride = 1 + x.shape[1] x.flat[ slice(0,None,stride) ] = v return x
Ooh that is clever. This should really go on the Cookbook page.
Perhaps more generally, for in place modification:
def setdiag(m, d): assert(len(x.shape) == 2) stride = 1 + x.shape[1] m.flat[slice(0,None,stride)] = d
Yeah, that's awesome. I have so many for-loops littering my code for setting diagonals. Here's an nd-version: def setdiag(a, d): assert(all([s == len(d) for s in a.shape])) stride = 1+sum(cumprod(a.shape[:-1])) a.flat[::stride] = d Cheers, Anand
On Wed, Jan 30, 2008 at 10:28 PM, Anand Patil<anand.prabhakar.patil@gmail.com> wrote:
Yeah, that's awesome. I have so many for-loops littering my code for setting diagonals. Here's an nd-version:
def setdiag(a, d): assert(all([s == len(d) for s in a.shape])) stride = 1+sum(cumprod(a.shape[:-1])) a.flat[::stride] = d
If someone feels like reviewing this ticket: http://projects.scipy.org/numpy/attachment/ticket/1132/numpy-index-funcs.dif... it has the above and a few more utilities, with docs and tests. Cheers, f
2009/7/1 Fernando Perez <fperez.net@gmail.com>:
If someone feels like reviewing this ticket:
http://projects.scipy.org/numpy/attachment/ticket/1132/numpy-index-funcs.dif...
it has the above and a few more utilities, with docs and tests.
Thanks! Reviewed and applied in r7101 through r7103. Cheers Stéfan
Hey Stefan, 2009/7/4 Stéfan van der Walt <stefan@sun.ac.za>:
2009/7/1 Fernando Perez <fperez.net@gmail.com>:
If someone feels like reviewing this ticket:
http://projects.scipy.org/numpy/attachment/ticket/1132/numpy-index-funcs.dif...
it has the above and a few more utilities, with docs and tests.
Thanks! Reviewed and applied in r7101 through r7103.
thanks for working on that ticket! BTW, looking at numpy/lib/index_tricks.py line 740: step = cumprod((1,)+a.shape[:-1]).sum() might be better written as suggested by Anand: step = 1+(cumprod(a.shape[:-1])).sum() which replaces a tuple concatenation by a simple numerical addition. Order-epsilon nit though. Thanks for the extra tests too! Best, f
On 30/01/2008, Carlos Scheidegger <cscheid@sci.utah.edu> wrote:
Say I have a matrix m whose diagonal values I want to replace with the contents of vector new_diagonal. Right now, I'm using
m -= diag(m) m += diag(new_diagonal)
If m is n by n, how about this? m[range(n),range(n)]=new_diagonal Slightly less efficient than using strides as suggested elsewhere if the array is contiguous in memory, but it'll work no matter how m is laid out and it's clear. It also generalizes, of course. (Incidentally, if it seems like this should set all elements of m, this is what I would have thought too, but it's not how numpy works. For that use the N._ix function.) Anne
To make up for my earlier response, I have coded a small timer routine that generates a 1000x1000 uniformly random matrix with values between 0 and 1 (using cvxopt), then tests the code suggested 10 000 times each with randomly generated diagonal vectors. I have a slow system, but the averaged results are: 0.8 seconds to run David's striding code (in place modification) 6 seconds to run Anne's code 18 seconds to run a for-loop This is an impressive speed up. For smaller matrices (100x100) there is still one order of magnitude between for-next loops and the suggested setdiag. Code available upon request (requires cvxopt). Cheers, Mclean
Just curious--is it safe to use the assert statement for anything beyond debugging in case someone actually runs with optimization? Cheers, William On Jan 31, 2008 3:00 AM, Mclean Edwards <mcleane@math.ubc.ca> wrote:
To make up for my earlier response, I have coded a small timer routine that generates a 1000x1000 uniformly random matrix with values between 0 and 1 (using cvxopt), then tests the code suggested 10 000 times each with randomly generated diagonal vectors.
I have a slow system, but the averaged results are:
0.8 seconds to run David's striding code (in place modification)
6 seconds to run Anne's code
18 seconds to run a for-loop
This is an impressive speed up.
For smaller matrices (100x100) there is still one order of magnitude between for-next loops and the suggested setdiag.
Code available upon request (requires cvxopt).
Cheers,
Mclean _______________________________________________ SciPy-user mailing list SciPy-user@scipy.org http://projects.scipy.org/mailman/listinfo/scipy-user
On 31/01/2008, william ratcliff <william.ratcliff@gmail.com> wrote:
Just curious--is it safe to use the assert statement for anything beyond debugging in case someone actually runs with optimization?
If debugging is turned off, then the condition in assert is not checked. So no. More importantly, AssertionErrors *mean* that there is a bug in your program. Not that it got invlaid input, not that it's out of some resource, they mean your program has a bug. So it's okay to turn them off if you're in a hurry, because a properly-working program never signals AssertionError no matter what you feed it. And really, how hard is it to replace assert c with if not c: raise ValueError ? Anne
On Sat, May 23, 2009 at 16:02, Alan G Isaac <aisaac@american.edu> wrote:
On 1/31/2008 1:37 AM Anne Archibald apparently wrote:
m[range(n),range(n)]=new_diagonal
Will that work with range objects (in Python 3)?
No. The automatic conversion to arrays does not consume iterators (nor will it when we port to Python 3). -- Robert Kern "I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth." -- Umberto Eco
On Sat, May 23, 2009 at 16:02, Alan G Isaac wrote:
On 1/31/2008 1:37 AM Anne Archibald apparently wrote:
m[range(n),range(n)]=new_diagonal Will that work with range objects (in Python 3)?
On 5/23/2009 5:05 PM Robert Kern apparently wrote:
No. The automatic conversion to arrays does not consume iterators (nor will it when we port to Python 3).
Sure, but range objects are not iterators. They are "almost" sequences. Python 3.0 (r30:67507, Dec 3 2008, 20:14:27) [MSC v.1500 32 bit (Intel)] on win32 Type "help", "copyright", "credits" or "license" for more information. >>> r = range(10) >>> next(r) Traceback (most recent call last): File "<stdin>", line 1, in <module> TypeError: range object is not an iterator >>> 5 in r True >>> list(r) [0, 1, 2, 3, 4, 5, 6, 7, 8, 9] >>> list(r) [0, 1, 2, 3, 4, 5, 6, 7, 8, 9] Cheers, Alan Isaac
On Sat, May 23, 2009 at 16:26, Alan G Isaac <aisaac@american.edu> wrote:
On Sat, May 23, 2009 at 16:02, Alan G Isaac wrote:
On 1/31/2008 1:37 AM Anne Archibald apparently wrote:
m[range(n),range(n)]=new_diagonal Will that work with range objects (in Python 3)?
On 5/23/2009 5:05 PM Robert Kern apparently wrote:
No. The automatic conversion to arrays does not consume iterators (nor will it when we port to Python 3).
Sure, but range objects are not iterators. They are "almost" sequences.
The answer is still no. Perhaps someone will write special support for that type when we do the Python 3 port, but there's nothing in numpy that would make it work automatically. For example, xrange() does not work as an index with the current numpy. -- Robert Kern "I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth." -- Umberto Eco
Robert Kern wrote:
On Sat, May 23, 2009 at 16:26, Alan G Isaac <aisaac@american.edu> wrote:
On Sat, May 23, 2009 at 16:02, Alan G Isaac wrote:
On 1/31/2008 1:37 AM Anne Archibald apparently wrote:
m[range(n),range(n)]=new_diagonal Will that work with range objects (in Python 3)? On 5/23/2009 5:05 PM Robert Kern apparently wrote: No. The automatic conversion to arrays does not consume iterators (nor will it when we port to Python 3). Sure, but range objects are not iterators. They are "almost" sequences.
The answer is still no. Perhaps someone will write special support for that type when we do the Python 3 port, but there's nothing in numpy that would make it work automatically. For example, xrange() does not work as an index with the current numpy.
Well, ranges are more capable than you think in Python 3: v = range(25) print (v[3], v[0], v[22], v) prints: 3 0 22 range(0, 25) -Scott
On Wed, Jul 1, 2009 at 12:21, Scott David Daniels<Scott.Daniels@acm.org> wrote:
Robert Kern wrote:
On Sat, May 23, 2009 at 16:26, Alan G Isaac <aisaac@american.edu> wrote:
On Sat, May 23, 2009 at 16:02, Alan G Isaac wrote:
On 1/31/2008 1:37 AM Anne Archibald apparently wrote:
m[range(n),range(n)]=new_diagonal Will that work with range objects (in Python 3)? On 5/23/2009 5:05 PM Robert Kern apparently wrote: No. The automatic conversion to arrays does not consume iterators (nor will it when we port to Python 3). Sure, but range objects are not iterators. They are "almost" sequences.
The answer is still no. Perhaps someone will write special support for that type when we do the Python 3 port, but there's nothing in numpy that would make it work automatically. For example, xrange() does not work as an index with the current numpy.
Well, ranges are more capable than you think in Python 3: v = range(25) print (v[3], v[0], v[22], v) prints: 3 0 22 range(0, 25)
No, they are exactly as capable as I think, i.e. as capable as xrange() is in Python 2: In [10]: v = xrange(25) In [11]: print v[3], v[0], v[22], v 3 0 22 xrange(25) Quite simply, numpy does not support arbitrary sequences and sequence-like objects as indices. If the eventual numpy port to Python 3 supports range() objects as indices, it will be because someone will have written special code for it. -- Robert Kern "I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth." -- Umberto Eco
participants (11)
-
Alan G Isaac -
Anand Patil -
Anne Archibald -
Carlos Scheidegger -
David Warde-Farley -
Fernando Perez -
Mclean Edwards -
Robert Kern -
Scott David Daniels -
Stéfan van der Walt -
william ratcliff