Наши проекты:
Журнал · Discuz!ML · Wiki · DRKB · Помощь проекту |
||
ПРАВИЛА | FAQ | Помощь | Поиск | Участники | Календарь | Избранное | RSS |
[18.119.110.116] |
|
Страницы: (2) [1] 2 все ( Перейти к последнему сообщению ) |
Сообщ.
#1
,
|
|
|
V etoj programmme modul sg1 v kotoroj funkcia dot
__author__ = 'kardashevskiy' import numpy as np from scipy import sparse from scipy.sparse import linalg from sg1 import sg1 import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D import time def f(x, y, r): return 2*c*(y-y**2+r*(x**2-x)) def poisson(n, h, r): V = [] I = [] J = [] b = [] for j in range(n): for i in range(n): k = i+n*j b.append(h**2*f(i*h+h, j*h+h, r)) V.append(2.-2*r) I.append(k) J.append(k) if i>0: V.append(-1.) I.append(k) J.append(k-1) if i<n-1: V.append(-1.) I.append(k) J.append(k+1) if j>0: V.append(r) I.append(k) J.append(k-n) if j<n-1: V.append(r) I.append(k) J.append(k+n) return sparse.coo_matrix((V, (I,J))), b N =100 h = 1./N c = 4. r = 0.4 tst = time.clock() A, b = poisson(N-1, h, r) B = A.tocsr() v, info = sg1(B, b,tol=1.e-08) dt = time.clock() - tst if iter == 0: print "N = %i, time solution = %1.3e" % (N, dt) fig = plt.figure() ax = Axes3D(fig) x = np.linspace(h/2, 1.-h/2, N-1) y = np.linspace(h/2, 1.-h/2, N-1) X,Y = np.meshgrid(x, y) Z = np.reshape(v, (N-1, N-1)) ax.plot_surface(X, Y, Z, rstride=2, cstride=2, cmap='spectral') plt.show() Vot modul __author__ = 'kardashevskiy' import numpy as np from scipy import sparse def sg1(B,b,tol=1.0e-8): """ Solve the linear system Ax=b by conjugate gradient method """ n=len(b) x=np.zeros((n),'float') r=np.copy(b) for i in range(n): r[i]=np.dot(B[i,0:n],x[0:n])-b[i] s=np.copy(r) As=np.zeros((n),'float') for k in range(n): for i in range(n): As[i]=np.dot(B[i,0:n],s[0:n]) alpha=np.dot(r,r)/np.dot(s,As) x=x-alpha*s for i in range(n): r[i]=np.dot(B[i,0:n],x[0:n])-b[i] if np.dot(r,r)<tol**2: break else: beta=-np.dot(r,As)/np.dot(s,As) s=r+beta*s return x,k A vot shto soobchaet /usr/bin/python /home/kardashevskiy/PycharmProjects/tests/qwe1.py Traceback (most recent call last): File "/home/kardashevskiy/PycharmProjects/tests/qwe1.py", line 47, in <module> v, info = sg1(B, b,tol=1.e-08) File "/home/kardashevskiy/PycharmProjects/tests/sg1.py", line 13, in sg1 r[i]=np.dot(B[i,0:n],x[0:n])-b[i] File "/usr/lib/python2.7/dist-packages/scipy/sparse/compressed.py", line 193, in __sub__ raise NotImplementedError('adding a scalar to a sparse ' NotImplementedError: adding a scalar to a sparse matrix is not supported Process finished with exit code 1 |
Сообщ.
#2
,
|
|
|
Навряд ли это от умножения выдаётся. Действительно вычитание скаляра не поддерживаемо разреженными матрицами, а вы b вычитаете. Надо или как-то указать, что результат будет обычная матрица, или менять сам метод.
dot тут не при чём. |
Сообщ.
#3
,
|
|
|
Что-то я прогоняю или не понимаю: там вообще строка из B на вектор множится, а потом число вычитается. В итоге всё в вектор уходит. Да, он не разрежен, но ошибку то пишет про матрицу, а не про вектор!
|
Сообщ.
#4
,
|
|
|
Я не знаю насколько правильно но я думаю так: В основной программе матрица V формируется в формате return sparse.coo_matrix((V, (I,J))), b то есть хранятся лишь координаты и значения ненулевых элементов матрицы. А dot ждет обычную запись. Я искал в этом направлении но не нашел
|
Сообщ.
#5
,
|
|
|
А просто попробуйте убрать вычитание b и увидите в чём проблема.
|
Сообщ.
#6
,
|
|
|
Я пробовал убирать это вычитание не получалось
Добавлено /usr/bin/python /home/kardashevskiy/PycharmProjects/tests/qwe1.py Traceback (most recent call last): File "/home/kardashevskiy/PycharmProjects/tests/qwe1.py", line 47, in <module> v, info = sg1(B, b,tol=1.e-08) File "/home/kardashevskiy/PycharmProjects/tests/sg1.py", line 13, in sg1 r[i]=np.dot(B[i,0:n],x[0:n]) ValueError: setting an array element with a sequence. Process finished with exit code 1 Добавлено А может подскажете как переформатировать с sparce cool в обычный формат чтобы dot заработал мне бы не хотелось менять основную программу |
Сообщ.
#7
,
|
|
|
У вас одинаковая строка в двух местах, а ошибку пишет про одно место:
r[i]=np.dot(B[i,0:n],x[0:n])-b[i] Попробуйте в одном месте написать так и посмотрим, в чём проблема: r[i]=-b[i] + np.dot(B[i,0:n],x[0:n]) |
Сообщ.
#8
,
|
|
|
Если так
__author__ = 'kardashevskiy' import numpy as np from scipy import sparse def sg1(B,b,tol=1.0e-8): """ Solve the linear system Ax=b by conjugate gradient method """ n=len(b) x=np.zeros((n),'float') r=np.copy(b) for i in range(n): r[i]=-b[i]+np.dot(B[i,0:n],x[0:n]) s=np.copy(r) As=np.zeros((n),'float') for k in range(n): for i in range(n): As[i]=np.dot(B[i,0:n],s[0:n]) alpha=np.dot(r,r)/np.dot(s,As) x=x-alpha*s for i in range(n): r[i]=np.dot(B[i,0:n],x[0:n]) if np.dot(r,r)<tol**2: break else: beta=-np.dot(r,As)/np.dot(s,As) s=r+beta*s return x,k Traceback (most recent call last): File "/home/kardashevskiy/PycharmProjects/tests/qwe1.py", line 47, in <module> v, info = sg1(B, b,tol=1.e-08) File "/home/kardashevskiy/PycharmProjects/tests/sg1.py", line 13, in sg1 r[i]=-b[i]+np.dot(B[i,0:n],x[0:n]) File "/usr/lib/python2.7/dist-packages/scipy/sparse/compressed.py", line 187, in __radd__ return self.__add__(other) File "/usr/lib/python2.7/dist-packages/scipy/sparse/compressed.py", line 173, in __add__ raise NotImplementedError('adding a scalar to a CSC or CSR ' NotImplementedError: adding a scalar to a CSC or CSR matrix is not supported /usr/bin/python /home/kardashevskiy/PycharmProjects/tests/qwe1.py Traceback (most recent call last): File "/home/kardashevskiy/PycharmProjects/tests/qwe1.py", line 47, in <module> v, info = sg1(B, b,tol=1.e-08) File "/home/kardashevskiy/PycharmProjects/tests/sg1.py", line 13, in sg1 r[i]=np.dot(B[i,0:n],x[0:n])-b[i] File "/usr/lib/python2.7/dist-packages/scipy/sparse/compressed.py", line 193, in __sub__ raise NotImplementedError('adding a scalar to a sparse ' NotImplementedError: adding a scalar to a sparse matrix is not supported Process finished with exit code 1 Добавлено если убрать совсем строку и т д то второй дот все равно находит |
Сообщ.
#9
,
|
|
|
Это программа прекрасно работает с библиотечным модулем linalg.cg который одинаково вычисляет с модулем с которым мы мучаемся в другой тестовой программе. То есть библиотечный модуль одинаково хорошо работает с любыми матрицами. Наш с вами модуль хорошо вычисляет в другой программе. Исходный код библиотечного модуля я не сумел найти. И код нашего модуля для меня более понятен и я хотел используя этот код попробовать похожие методы модули так вот приведу код который работает.
__author__ = 'kardashevskiy' #Hyperbolic equation import numpy as np from scipy import sparse from scipy.sparse import linalg import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D import time def f(x, y, r): return 2*c*(y-y**2+r*(x**2-x)) def poisson(n, h, r): V = [] I = [] J = [] b = [] for j in range(n): for i in range(n): k = i+n*j b.append(h**2*f(i*h+h, j*h+h, r)) V.append(2.-2*r) I.append(k) J.append(k) if i>0: V.append(-1.) I.append(k) J.append(k-1) if i<n-1: V.append(-1.) I.append(k) J.append(k+1) if j>0: V.append(r) I.append(k) J.append(k-n) if j<n-1: V.append(r) I.append(k) J.append(k+n) return sparse.coo_matrix((V, (I,J))), b N =100 h = 1./N c = 4. r = 0.4 tst = time.clock() A, b = poisson(N-1, h, r) B = A.tocsr() v, info = linalg.cg(B, b, tol=1.e-08) dt = time.clock() - tst if iter == 0: print "N = %i, time solution = %1.3e" % (N, dt) fig = plt.figure() ax = Axes3D(fig) x = np.linspace(h/2, 1.-h/2, N-1) y = np.linspace(h/2, 1.-h/2, N-1) X,Y = np.meshgrid(x, y) Z = np.reshape(v, (N-1, N-1)) ax.plot_surface(X, Y, Z, rstride=2, cstride=2, cmap='spectral') plt.show() |
Сообщ.
#10
,
|
|
|
Цитата Urung @ А давайте так попробуем: вместо одной итерации:если убрать совсем строку и т д то второй дот все равно находит for i in range(n): r[i]=np.dot(B[i,0:n],x[0:n])-b[i] for i in range(n): r[i] = np.dot(B[i,0:n],x[0:n]) r[i] = r[i]-b[i] Вот тогда и посмотрим на что будет ругаться, на вычитание или на умножение. Добавлено Ай, тьфу, вы ж в №6 написали, что не работает так. Виноват. Думаю. Добавлено Хм... смотрите, интересный момент: 1.Тут пишет ошибку: Цитата Urung @ File "/home/kardashevskiy/PycharmProjects/tests/sg1.py", line 13, in sg1 r[i]=np.dot(B[i,0:n],x[0:n]) ValueError: setting an array element with a sequence. 2.а тут норма: As[i]=np.dot(B[i,0:n],s[0:n]) |
Сообщ.
#11
,
|
|
|
Забавно, но по этому коду:
Цитата Urung @ можно сделать вывод, что r = -b, т.к. x - нулевой. x=np.zeros((n),'float') r=np.copy(b) for i in range(n): r[i]=np.dot(B[i,0:n],x[0:n])-b[i] |
Сообщ.
#12
,
|
|
|
СлавянСпасибо за участие но тут кажется тупик. Но вот может быть посоветуете найти исходный код библиотечного ьодуля linalg.cg я пытался искал но не смог найти я попробовал бы преобразовать тот модуль наверняка для вас это не проблема наверное вы вовсю используете модифицируете и новые делаете. А мы чайники даже не можем найти.
|
Сообщ.
#13
,
|
|
|
Навряд ли "посоветую". Я ни в питоне ни в Руби "не шарю", а вам советую создать тему в разделе С/С++ Алгоритмы про linalg, а там много умного народа бывает, они явно укажут. Даже если найдёте не питоновский, то переделать на него будет не шибко сложно. Так что "рвите когти" туда!
|
Сообщ.
#14
,
|
|
|
Спасибо я вам очень признателен за то что приняли очень живое участие. Спасибо попробую надеюсь не побьют.
|
Сообщ.
#15
,
|
|
|
Urung, ты алгоритм с какого-то другого языка переписывал? Или может в какой-то книге, где он через вычисления с матрицами матрицы был записан?
Если так, то надо учитывать, что например в MathCAD'е или MATLAB'е матрица размером 1х1 может быть использована там, где нужен скаляр (отдельное число). В Python'овских библиотеках numpy и scipy не так - матрица остаётся матрицей, даже если она вообще не содержит элементов или содержит единственный элемент. В частности чтобы получить из такой одноэлементной матрицы её единственный элемент в пакете numpy есть функция asscalar. Можешь попробовать написать r[i] = np.asscalar(np.dot(B[i,0:n],x[0:n])) - b[i] |