pysparse-users

[Pysparse-users] PysparseSpDiagsMatrix usage from spdiags in matlab

[zipeppe <zi...@gm...>]

Dear all,

I don't know it this is the right place to ask, but unfortunately the
official Python-forum doesn't help.

I am quite new to Python and I am unable to correctly translate the
following MATLAB script in Python code, specially because of the
*spdiags*proprietary function:

function bd=BLC(b,A,B,s)
   L=length(b);
   Bs=-B/s;  As=-A/s;
   bd=ones(L,1)*median(b);bd0=b;nm=norm(bd-bd0);nm0=RealMax;
   M0=-ones(L,1)/As;
   e=ones(L,1);
   D0=spdiags([1*e -4*e 6*e -4*e 1*e],-2:2,L,L);
   D0(1,1)=2; D0(L,L)=2;
   D0(2,1)=-4; D0(1,2)=-4; D0(L,L-1)=-4; D0(L-1,L)=-4;
   D0(2,2)=10;D0(L-1,L-1)=10;
   I=0;
   while nm>10 & I<30
       %& nm<nm0;
       I=I+1;
       M=M0;D=D0;bd0=bd;nm0=nm;
       for i=1:L
           if bd(i)>b(i)
               M(i)=M(i)+2*Bs*b(i)/As;
               D(i,i)=D(i,i)+2*Bs/As;
           end
       end
       bd=D\M;
       nm=norm(bd0-bd);
   end

The algorithm receives *b* as (Nx1 data vector), *A* (1x1), *B* (1x1) and *s
* (1x1) as local parameters, and creates *bd* (Nx1data vector) as output.

It looks like the PySparse package is the best promising, despite that I am
not able to use *PysparseSpDiagsMatrix* function in any way... I got only
errors like:

*TypeError: unsubscriptable object  *
- or -
*IndexError: Not as many row indices as values*

Please, help me to understand what is going on in MATLAB and how to
translate this operation using the PySparse module:

Let have a B matrix as follow:

B =
   1    6   11
   2    7   12
   3    8   13
   4    9   14
   5   10   15

The spdiags function works as follows:

A = full (spdiags (B, [-2 0 2], 5, 5))

 Matrix B                      Matrix A

1    6   11                 6    0   13    0    0
2    7   12                 0    7    0   14    0
3    8   13  == spdiags =>  1    0    8    0   15
4    9   14                 0    2    0    9    0
5   10   15                 0    0    3    0   10

A(3,1), A(4,2), and A(5,3) are taken from the upper part of B(:,1)

A(1,3), A(2,4), and A(3,5) are taken from the lower part of B(:,3).
How can I implement the same operation in python using *
PysparseSpDiagsMatrix* ?



Thanks for your help,
ZPP

-- 

-------------------------------------------------------------------
Coltivate Linux che Windows si pianta da solo!
-------------------------------------------------------------------

[Pysparse-users] PysparseSpDiagsMatrix usage from spdiags in matlab

[zipeppe <zi...@gm...>]

Dear all,

I don't know it this is the right place to ask, but unfortunately the
official Python-forum doesn't help.

I am quite new to Python and I am unable to correctly translate the
following MATLAB script in Python code, specially because of the
*spdiags*proprietary function:

function bd=BLC(b,A,B,s)
   L=length(b);
   Bs=-B/s;  As=-A/s;
   bd=ones(L,1)*median(b);bd0=b;nm=norm(bd-bd0);nm0=RealMax;
   M0=-ones(L,1)/As;
   e=ones(L,1);
   D0=spdiags([1*e -4*e 6*e -4*e 1*e],-2:2,L,L);
   D0(1,1)=2; D0(L,L)=2;
   D0(2,1)=-4; D0(1,2)=-4; D0(L,L-1)=-4; D0(L-1,L)=-4;
   D0(2,2)=10;D0(L-1,L-1)=10;
   I=0;
   while nm>10 & I<30
       %& nm<nm0;
       I=I+1;
       M=M0;D=D0;bd0=bd;nm0=nm;
       for i=1:L
           if bd(i)>b(i)
               M(i)=M(i)+2*Bs*b(i)/As;
               D(i,i)=D(i,i)+2*Bs/As;
           end
       end
       bd=D\M;
       nm=norm(bd0-bd);
   end

The algorithm receives *b* as (Nx1 data vector), *A* (1x1), *B* (1x1) and *s
* (1x1) as local parameters, and creates *bd* (Nx1data vector) as output.

It looks like the PySparse package is the best promising, despite that I am
not able to use *PysparseSpDiagsMatrix* function in any way... I got only
errors like:

*TypeError: unsubscriptable object  *
- or -
*IndexError: Not as many row indices as values*

Please, help me to understand what is going on in MATLAB and how to
translate this operation using the PySparse module:

Let have a B matrix as follow:

B =
   1    6   11
   2    7   12
   3    8   13
   4    9   14
   5   10   15

The spdiags MATLAB function works as follows:

A = full (spdiags (B, [-2 0 2], 5, 5))

 Matrix B                      Matrix A

1    6   11                 6    0   13    0    0
2    7   12                 0    7    0   14    0
3    8   13  == spdiags =>  1    0    8    0   15
4    9   14                 0    2    0    9    0
5   10   15                 0    0    3    0   10

A(3,1), A(4,2), and A(5,3) are taken from the upper part of B(:,1)

A(1,3), A(2,4), and A(3,5) are taken from the lower part of B(:,3).
How can I implement the same operation in python using *
PysparseSpDiagsMatrix* ?



Thanks for your help,
ZPP


-- 

-------------------------------------------------------------------
Coltivate Linux che Windows si pianta da solo!
-------------------------------------------------------------------

Re: [Pysparse-users] PysparseSpDiagsMatrix usage from spdiags in matlab

[Dominique O. <dom...@gm...>]

On Thu, Nov 5, 2009 at 4:00 AM, zipeppe <zi...@gm...> wrote:

> Dear all,
>
> I don't know it this is the right place to ask, but unfortunately the
> official Python-forum doesn't help.
>

This is the right place for Pysparse questions.


> I am quite new to Python and I am unable to correctly translate the
> following MATLAB script in Python code, specially because of the *spdiags*proprietary function:
>
> function bd=BLC(b,A,B,s)
>    L=length(b);
>    Bs=-B/s;  As=-A/s;
>    bd=ones(L,1)*median(b);bd0=b;nm=norm(bd-bd0);nm0=RealMax;
>    M0=-ones(L,1)/As;
>    e=ones(L,1);
>    D0=spdiags([1*e -4*e 6*e -4*e 1*e],-2:2,L,L);
>    D0(1,1)=2; D0(L,L)=2;
>    D0(2,1)=-4; D0(1,2)=-4; D0(L,L-1)=-4; D0(L-1,L)=-4;
>    D0(2,2)=10;D0(L-1,L-1)=10;
>    I=0;
>    while nm>10 & I<30
>        %& nm<nm0;
>        I=I+1;
>        M=M0;D=D0;bd0=bd;nm0=nm;
>        for i=1:L
>            if bd(i)>b(i)
>                M(i)=M(i)+2*Bs*b(i)/As;
>                D(i,i)=D(i,i)+2*Bs/As;
>            end
>        end
>        bd=D\M;
>        nm=norm(bd0-bd);
>    end
>
> The algorithm receives *b* as (Nx1 data vector), *A* (1x1), *B* (1x1) and
> *s* (1x1) as local parameters, and creates *bd* (Nx1data vector) as
> output.
>

The Matlab command that you highlighted creates an L x L matrix with the
vectors e, -4*e, 6*e, -4*e and e on and around the main diagonal. The
position argument is -2:2 which means that the first vector (e) goes on the
sub-sub-diagonal, the second vector (-4e) goes on the sub-diagonal, 6e goes
on the main diagonal, -4e goes one above the main diagonal and e goes two
above the main diagonal like so.

Here is a code fragment that performs the same operation with Pysparse:

In [1]: from pysparse.pysparseMatrix import PysparseSpDiagsMatrix as spdiags
In [2]: import numpy as np
In [3]: L = 5  # The length of b in your Matlab code
In [4]: e = np.ones(L)
In [5]: D0 = spdiags(L, (e, -4*e, 6*e, -4*e, e), (-2,-1,0,1,2))
In [6]: print D0
 6.000000  -4.000000   1.000000      ---        ---
-4.000000   6.000000  -4.000000   1.000000      ---
 1.000000  -4.000000   6.000000  -4.000000   1.000000
    ---     1.000000  -4.000000   6.000000  -4.000000
    ---        ---     1.000000  -4.000000   6.000000

The first argument (L) specifies the matrix size, the second is the set of
vectors that should be laid out on the diagonals and the third gives the
indices of the diagonals that are concerned.

For a more Matlab-like notation, the last argument can also be written
np.r_[-2:2] (see the Numpy function r_ ... note the underscore).

I hope this helps.

-- 
Dominique