Erstellung von Index-Arrays für csc_matrix

[Shaint]

Hallo,
ich versuche momentan ein Programm das ich in Octave geschrieben habe in Python zu übersetzten.
Das Programm berechnet am Ende ein Sparse-Eigenwertproblem. Ich scheitere allerdings bereits beim Erstellen meiner csc_matrix.
Ich möchte mir ein Daten-Array, ein Column _Index- und ein Row_Index-Array erzeugen und dies nutzen um meine Martix zu erzeugen.
Beispielhaft möchte ich eine 5x5 große Matrix erzeugen die mit null initialisiert ist und einem Rand r der Breite r=2 den Wert 5000 zuweist.
Das klappt auch soweit für den oberen und den unteren Rand, allerdings sobald ich mehr als zwei Anweisungen der Form: [np.arange(r)]*ny (ny ist die Spaltenanzahl) in die Anweisung np.array() schreibe bekomme ich eine Fehlermeldung.
Meine Frage ist nun, wie erzeuge ich ein Array (möglichst in einer Zeile) das mehrer Anweisungen der Form [np.arange(r)]*ny verkraftet.
Damit es deutlicher wird hier der Octave-Code der 2 Matrizen die ich erzeugen möchte:
Matrix 1:

Code: Alles auswählen

nx = 4
	ny = 5
n=nx*ny;
r=1;
sigma_Cu = 5e7;
sigma_mat = zeros(nx,ny);
sigma_mat(1:r,:) = sigma_Cu; 
sigma_mat(end:-1:end-(r-1),:) = sigma_Cu; 
sigma_mat(:,1:r) = sigma_Cu; 
sigma_mat(:,end:-1:end-(r-1)) = sigma_Cu;

Diese Matrix ist nur am Rand ungleich Null.
Für den oberen Rand funktioniert das auch noch ganz gut in Python:

Code: Alles auswählen

x = 5000
times = ny*r*2+(nx-2*r)*2 
data_vec = np.array([x]*ny*r)
A_row = np.array([np.arange(r)]*ny).flatten()
A_col = np.array([np.arange(ny)]*r).flatten()
A = scipy.sparse.csc_matrix((data_vec,(A_row,A_col)),shape=[nx,ny])

Sobald ich nun versuche bpsw. die Index-Arrays für den linken Rand anzuhängen an A_row und A_col kommt es zu einer Reihe Fehlermeldungen.

Matrix 2:

Code: Alles auswählen

diagvec = ones(1,n);
diagvec2 = ones(1,n);
diagvec2(1:nx:end) = 2;
diagvec3 = ones(1,n);
diagvec3(nx:nx:end) = 0;

A = sparse([1:n],[1:n],diagvec,n,n);
A = A + sparse([1:nx],[1:nx],diagvec(1:nx),n,n);
A = A + sparse([nx+1:n],[1:n-nx],diagvec(1:n-nx),n,n);

(sparse(Zeilenindices,Spaltenindices,Daten,Zeilenanzahl,Spaltenanzahl))
Bei der Erzeugung dieser Matrix in Python scheitert es ebenfalls daran die notwendigen Arrays zu erzeugen.
Ich dachte nun, dass ich um diese Matrix zu realisieren nur die Spalten-Indices und die Zeilen-Indices in jeweils ein Array zusammenfassen muss und die Datenarrays (diavec-diavec3) ebenfalls in ein Array zusammen fasse.
Ich wäre euch sehr dankbar, wenn ihr mir erklären könntet wie ich das bewerkstellige.
Ich hoffe ich hab mich verständlich ausgedrückt


Grüße
Shaint

[EyDu]

Hallo und willkommen im Forum!

Ich habe mir den Quelltext jetzt nicht angeschaut, aber wie sieht denn die Fehlermeldung, inklusive des gesamten Tracebacks, aus? Die geben häufig einen Hinweis auf den Fehler

[Shaint]

Ja ähm, das war etwas dümmlich den Fehler-Code nicht mit zu teilen

Traceback (most recent call last):
File "H:\.temp.py", line 55, in <module>
A_row = np.array([np.arange(r)]*ny,[np.arange(nx-r,nx)]*ny).flatten()
TypeError: data type not understood

(np: import numpay as np)
Die einzelnen Arange-Arrays werden wie gewünscht erzeugt, allerdings eben nicht in einem gemeinsamen array verknüpft. Da noch drei solche arange-Arrays kommen und mehrere derartige Matrizen möchte ich das ungern immer in einer for-Schleife mit np.concatenate realisieren.

[BlackJack]

@Shaint: Das Komma separiert verschiedene Argumente und `np.array()` möchte als zweites (optionales) Argument den Datentyp haben, den das erstellte Array haben soll. Und Du übergibst als Typ ein Array, das ist aber kein Datentyp mit dem `np.array()` etwas anfangen kann.

Die `concatenate()`-Funktion hast Du ja schon gefunden. Du brauchst da doch keine Schleife für. Die Funktion nimmt eine Sequenz mit allen Arrays entgegen, die zusammengefügt werden sollen.

[Shaint]

Danke für die Antwort.
Ich habe die Maitrix nun realisiert:

Code: Alles auswählen

sigma_Cu = 5e7
Cu_row_data = array([sigma_Cu] * (ny*r) )
Cu_col_data = array([sigma_Cu] * ((nx-2*r)*r) )

sigma_mat = sp.csc_matrix( (nx,ny))
row = array([arange(0,r)]*ny).flatten(1)
col = array([arange(0,ny)]*r).flatten()
sigma_mat = sigma_mat + sp.csc_matrix((Cu_row_data,(row,col)), shape=(nx,ny))
row = array([arange(nx-r,nx)]*ny).flatten(1)
col = array([arange(0,ny)]*r).flatten()
sigma_mat = sigma_mat + sp.csc_matrix((Cu_row_data,(row,col)), shape=(nx,ny))
row = array([arange(r,nx-r)]*r).flatten(1)
col = array([arange(0,r)]*(nx-2*r)).flatten()
sigma_mat = sigma_mat + sp.csc_matrix((Cu_col_data,(row,col)), shape=(nx,ny))
row = array([arange(r,nx-r)]*r).flatten(1)
col = array([arange(ny-r,ny)]*(nx-2*r)).flatten()
sigma_mat = sigma_mat + sp.csc_matrix((Cu_col_data,(row,col)), shape=(nx,ny))

Das Ergeniss ist auch zu meiner vollen Zufriedenheit.
Jedoch wie man sieht brauchte ich nun 16 Zeilen Code in Python um etwas zu realisieren, das ich in Octave in 5 Zeilen geschafft habe.
Die 2.Matrix muss ich noch erstellen, sowie viele weitere die alle in noch mehr Schritten entstehen (bis zu 5 Schritte). Gibt es vll. ein Tutorial, das sich mit der Erstellung von solchen Matrizen beschäftigt?
Mir gefällt die Funktinon "scipy.sparse.spdiags" recht gut, da sich mit ihr die Matrizen sehr einfach erstellen lassen, jedoch werden meine Matrizen sehr(!) groß und sind sehr dünn besetzt (sie besitzen eine dünne Bandstruktur). Am Ende soll ein Eigenwertproblem gelöst werden und der Löser (Arnoldi) benötigt dafür csc-Matrizen. Gibt es eine Möglichkeit Matrizen als csc zu erstellen, aber dennoch eine kompakte Schreibweise wie sparse.spdiags anzuwenden?

Ich bin euch für eure Hilfe sehr dankbar.

[BlackJack]

@Shaint: Die Dokumentation zu `scypi.sparse` legt ja nahe, dass man `lil_matrix` zum Aufbau verwendet und dann erst zum Rechnen in eine andere interne Darstellung umwandelt. Also zum Beispiel:

Code: Alles auswählen

from scipy.sparse import lil_matrix


def main():
    nx, ny = 6, 6
    r = 2
    sigma_Cu = 5e7

    sigma_mat = lil_matrix((nx, ny))
    sigma_mat[:r, :] = sigma_Cu
    sigma_mat[-r:, :] = sigma_Cu
    sigma_mat[:, :r] = sigma_Cu
    sigma_mat[:, -r:] = sigma_Cu
    sigma_mat = sigma_mat.tocsc()

    print sigma_mat.todense()


if __name__ == '__main__':
    main()

[Shaint]

Für die erste Matrix ist lil_matrix wirklich geschickter. Die anderen Matritzen gingen tatsächlich ganz einfach mit der standard csc_matrix. Wenn mann mit [a:b] richtig indiziert .
Habe mich da einfach dumm angestellt und hätte vll. noch öfter nachlesen sollen.
Vielen Dank für deine Hilfe und Geduld, nächstes mal stell ich mich etwas geschickter an und frag nicht gleich nach nur weils nach einigen Versuchen nicht klappt.

Nochmals vielen Dank!