Форум сайта python.su
Добрый вечер!
Пытаюсь решить мотодом наименьших квадратов систему двух матриц Ab=0.
pts2d = np.array([[1.0486, -0.3645], [-1.6851, -0.4004], [-0.9437, -0.42], [1.0682, 0.0699], [0.6077, -0.0771], [1.2543, -0.6454], [-0.2709, 0.8635], [-0.4571, -0.3645], [-0.7902, 0.0307], [0.7318, 0.6382], [-1.058, 0.3312], [0.3464, 0.3377], [0.3137, 0.1189], [-0.431, 0.0242], [-0.4799, 0.292], [0.6109, 0.083], [-0.4081, 0.292], [-0.1109, -0.2992], [0.5129, -0.0575], [0.1406, -0.4527]]) pts3d = np.array([[1.5706, -0.149, 0.2598], [-1.5282, 0.9695, 0.3802], [-0.6821, 1.2856, 0.4078], [0.4124, -1.0201, -0.0915], [1.2095, 0.2812, -0.128], [0.8819, -0.8481, 0.5255], [-0.9442, -1.1583, -0.3759], [0.0415, 1.3445, 0.324], [-0.7975, 0.3017, -0.0826], [-0.4329, -1.4151, -0.2774], [-1.1475, -0.0772, -0.2667], [-0.5149, -1.1784, -0.1401], [0.1993, -0.2854, -0.2114], [-0.432, 0.2143, -0.1053], [-0.7481, -0.384, -0.2408], [0.8078, -0.1196, -0.2631], [-0.7605, -0.5792, -0.1936], [0.3237, 0.797, 0.217], [1.3089, 0.5786, -0.1887], [1.2323, 1.4421, 0.4506]]) u = pts2d[:,0] v = pts2d[:,1] X = pts3d[:,0] Y = pts3d[:,1] Z = pts3d[:,2] A = np.zeros(40,12) for i in range[1:2:40]: for j in range ((i + 1) / 2): A([i:i + 1], [:]) = mcat([X(j), Y(j), Z(j) [ 1, 0, 0, 0, 0] - u(j)*X(j)-u(j)*Y(j)-u(j)*Z(j)-u(j), [0, 0, 0, 0] X(j), Y(j), Z(j), [1] - v(j)*X(j)-v(j)*Y(j)-v(j)*Z(j)-v(j)])
Офлайн
A = np.zeros(40,12) for i in range[1:2:40]: # перед этим for не нужен отступ for j in range ((i + 1) / 2): # отступ 4 пробела # в строке ниже отступ 8 пробелов A([i:i + 1], [:]) = mcat([X(j), Y(j), Z(j) [ 1, 0, 0, 0, 0] - u(j)*X(j)-u(j)*Y(j)-u(j)*Z(j)-u(j), [0, 0, 0, 0] X(j), Y(j), Z(j), [1] - v(j)*X(j)-v(j)*Y(j)-v(j)*Z(j)-v(j)])
Офлайн
Циклы в Python медленные, поэтому для решения систем линейных уравнений и прочих операций линейной алгебры, крайне желательно использовать numpy.
МНК решение системы
Ax=b
будет выглядеть (с использованием numpy)
import numpy as np x = np.linalg.pinv(A) * b
Офлайн
Спасибо, отступы исправила, но проблема была не в этом.
numpy использую, но для этой задачи не подходит, т.к. нужно прописать A вручную
pts2d = np.array([[1.0486, -0.3645], [-1.6851, -0.4004], [-0.9437, -0.42], [1.0682, 0.0699], [0.6077, -0.0771], [1.2543, -0.6454], [-0.2709, 0.8635], [-0.4571, -0.3645], [-0.7902, 0.0307], [0.7318, 0.6382], [-1.058, 0.3312], [0.3464, 0.3377], [0.3137, 0.1189], [-0.431, 0.0242], [-0.4799, 0.292], [0.6109, 0.083], [-0.4081, 0.292], [-0.1109, -0.2992], [0.5129, -0.0575], [0.1406, -0.4527]]) pts3d = np.array([[1.5706, -0.149, 0.2598], [-1.5282, 0.9695, 0.3802], [-0.6821, 1.2856, 0.4078], [0.4124, -1.0201, -0.0915], [1.2095, 0.2812, -0.128], [0.8819, -0.8481, 0.5255], [-0.9442, -1.1583, -0.3759], [0.0415, 1.3445, 0.324], [-0.7975, 0.3017, -0.0826], [-0.4329, -1.4151, -0.2774], [-1.1475, -0.0772, -0.2667], [-0.5149, -1.1784, -0.1401], [0.1993, -0.2854, -0.2114], [-0.432, 0.2143, -0.1053], [-0.7481, -0.384, -0.2408], [0.8078, -0.1196, -0.2631], [-0.7605, -0.5792, -0.1936], [0.3237, 0.797, 0.217], [1.3089, 0.5786, -0.1887], [1.2323, 1.4421, 0.4506]]) u = pts2d[:,0] v = pts2d[:,1] X = pts3d[:,0] Y = pts3d[:,1] Z = pts3d[:,2] A = np.zeros(40,12) for i in range(1,2,40): for j in range ((i + 1) / 2): A([i:i + 1], :) = np.concatenate(X(j), Y(j), Z(j) 1, 0, 0, 0, 0 - u(j)*X(j)-u(j)*Y(j)-u(j)*Z(j)-u(j) 0, 0, 0, 0 X(j), Y(j), Z(j), 1 - v(j)*X(j)-v(j)*Y(j)-v(j)*Z(j)-v(j)) U, D, V = np.linalg.svd(A) F = V[-1].reshape(3,3)
Отредактировано Rijdahl (Ноя. 8, 2016 09:34:46)
Офлайн
Rijdahl
Вы отступы не исправили, только еще хуже все сделали, просто скажите зачем Вы на питоне это пытаетесь сделать? Вы можете написать это на знакомом Вам языке?
потом:
>>> list(range(1,2,40)) [1]
A([i:i + 1], :)
Офлайн
Потому что передо мной поставлена именно такая задача = сделать в питоне До этого работала только в MatLab и ещё не хватает знания языка.
Исправила предыдущую ошибку в A:
A[i:i+1, :] = np.concatenate([ X[j], Y[j], Z[j], 1, 0, 0, 0, 0 - u[j]*X[j]-u[j]*Y[j]-u[j]*Z[j]-u[j], 0, 0, 0, 0, X[j], Y[j], Z[j], 1 - v[j]*X[j]-v[j]*Y[j]-v[j]*Z[j]-v[j] ])
Офлайн
напишите как бы Вы сделали в Матлабе
я только восстанавливаю свою математику
Офлайн
Первую функцию я уже сделала. Там действительно нужно было поменять A. Вот что получилось:
def createSystemForP(): pts2d = np.array([[1.0486, -0.3645], [-1.6851, -0.4004], [-0.9437, -0.42], [1.0682, 0.0699], [0.6077, -0.0771], [1.2543, -0.6454], [-0.2709, 0.8635], [-0.4571, -0.3645], [-0.7902, 0.0307], [0.7318, 0.6382], [-1.058, 0.3312], [0.3464, 0.3377], [0.3137, 0.1189], [-0.431, 0.0242], [-0.4799, 0.292], [0.6109, 0.083], [-0.4081, 0.292], [-0.1109, -0.2992], [0.5129, -0.0575], [0.1406, -0.4527]]) pts3d = np.array([[1.5706, -0.149, 0.2598], [-1.5282, 0.9695, 0.3802], [-0.6821, 1.2856, 0.4078], [0.4124, -1.0201, -0.0915], [1.2095, 0.2812, -0.128], [0.8819, -0.8481, 0.5255], [-0.9442, -1.1583, -0.3759], [0.0415, 1.3445, 0.324], [-0.7975, 0.3017, -0.0826], [-0.4329, -1.4151, -0.2774], [-1.1475, -0.0772, -0.2667], [-0.5149, -1.1784, -0.1401], [0.1993, -0.2854, -0.2114], [-0.432, 0.2143, -0.1053], [-0.7481, -0.384, -0.2408], [0.8078, -0.1196, -0.2631], [-0.7605, -0.5792, -0.1936], [0.3237, 0.797, 0.217], [1.3089, 0.5786, -0.1887], [1.2323, 1.4421, 0.4506]]) u = pts2d[:,0] v = pts2d[:,1] X = pts3d[:,0] Y = pts3d[:,1] Z = pts3d[:,2] A = np.zeros((40,12)) for i in range(0, 40, 2): for j in range (20): A[i, :] = np.array([ X[j], Y[j], Z[j], 1, 0, 0, 0, 0, - u[j]*X[j],-u[j]*Y[j],-u[j]*Z[j],-u[j]]) A[i+1, :] = np.array([ 0, 0, 0, 0, X[j], Y[j], Z[j], 1, - v[j]*X[j],-v[j]*Y[j],-v[j]*Z[j],-v[j]]) return A
[U D V] = svd(A); disp(V) l=size(V,1); m=V(l,:); disp(m) l=1; M=zeros(3,3); for i=1:3 for j=1:3 M(i,j) = m(l); l=l+1; end end
def solveForP(A): (U, D, V) = np.linalg.svd(A) print(V) l = len(V, 1) m = V[l, :] print(m) l = 1 P = np.zeros(3, 3) for i in range (1, 3, 1): for j in range (1, 3, 1): P[i, j]= m[i] l = l + 1 return P
Офлайн
%%matlab -i A -o M
[U D V] = svd(A);
disp(V)
l=size(V,1);
m=V(l,:);
disp(m)
l=1;
M=zeros(3,3);
for i=1:3
for j=1:3
M(i,j) = m(l);
l=l+1;
end
end
0.1507 0.5371 -0.0914 0.4104 0.1214 -0.4653 0.0360 -0.0180 -0.2094 0.1926 -0.4376 -0.0819
0.1764 0.6285 0.0976 -0.0856 0.0143 0.1157 0.0350 0.0430 0.0632 0.0120 0.6512 -0.3348
0.0551 0.1964 0.1370 -0.3198 -0.0397 0.4021 -0.0789 0.0336 0.1464 -0.2315 -0.5957 -0.4867
0.1223 0.4358 0.1163 -0.1921 -0.0018 0.2496 -0.0775 0.0211 0.0863 -0.0927 -0.1326 0.8014
-0.4853 0.1668 0.0194 -0.5878 -0.1291 -0.5563 0.2106 -0.0016 0.1350 -0.0347 -0.0264 0.0102
-0.5679 0.1952 -0.4003 0.4334 0.0392 0.1752 -0.0252 0.0252 0.3781 -0.3378 0.0242 0.0091
-0.1775 0.0610 0.4680 0.2856 -0.6740 0.0479 0.0245 -0.4532 -0.0213 0.0318 -0.0086 -0.0032
-0.3938 0.1354 -0.3210 -0.1782 -0.0785 0.3331 -0.1755 -0.0750 -0.6549 0.3300 0.0053 -0.0190
-0.2409 0.0000 0.3592 -0.0168 0.5603 -0.0658 -0.5283 -0.4101 0.1421 0.1580 0.0322 -0.0296
-0.2819 0.0000 0.4808 0.1701 -0.0409 -0.0284 -0.2184 0.7678 -0.1426 0.0174 0.0071 0.0000
-0.0881 0.0000 0.0056 0.0264 0.0092 0.2230 0.2881 0.1177 0.4615 0.7891 -0.0959 0.0000
-0.1955 0.0000 0.3269 0.0923 0.4348 0.1945 0.7065 -0.1140 -0.2814 -0.1655 -0.0153 0.0249
-0.1955 0.0000 0.3269 0.0923 0.4348 0.1945 0.7065 -0.1140 -0.2814 -0.1655 -0.0153 0.0249
print(matrix_format(M))
-0.1955 0.0000 0.3269
0.0923 0.4348 0.1945
0.7065 -0.1140 -0.2814
def row_format(row): return ' '.join('% .4f' % val for val in row) def matrix_format(mat): return '\n'.join(row_format(row) for row in mat) def solveForP(A): (U, D, V) = np.linalg.svd(A, full_matrices=1, compute_uv=1) print(matrix_format(V)) m = V[-1] print() print(row_format(m)) return m[:9].reshape((3, 3)) m = solveForP(A)
-0.1507 -0.1764 -0.0551 -0.1223 0.4853 0.5679 0.1775 0.3938 0.2409 0.2819 0.0881 0.1955
0.5371 0.6285 0.1964 0.4358 0.1668 0.1952 0.0610 0.1354 -0.0000 -0.0000 -0.0000 -0.0000
0.0597 0.4611 -0.3154 -0.6248 -0.1231 0.3748 -0.0189 -0.2886 0.0837 -0.2180 0.0154 0.0019
-0.0365 -0.1118 0.0519 0.1012 0.4392 0.2638 -0.2463 -0.5479 -0.5748 0.1218 0.0304 -0.0613
-0.0517 0.0261 -0.0790 -0.1607 0.0705 0.0230 0.3074 0.4573 -0.5335 -0.2310 -0.0710 -0.5590
0.0013 -0.0533 0.0409 0.0813 -0.5765 0.4364 -0.4645 0.2119 -0.1147 0.3812 -0.0565 -0.2095
0.0941 0.0179 -0.6710 0.1970 0.2275 -0.1776 -0.1556 -0.0713 0.2981 0.2555 0.0166 -0.4845
0.0268 -0.1093 -0.1447 0.2304 -0.3440 0.2121 0.7461 -0.3490 -0.0481 0.2619 -0.0115 -0.0240
-0.0232 -0.2888 -0.2055 0.4135 -0.0698 0.3825 -0.0984 -0.0211 0.1406 -0.7207 -0.0487 0.0041
-0.0892 -0.0528 0.5762 -0.0959 0.0939 0.1246 0.0544 -0.2483 0.4423 -0.0505 -0.0758 -0.5960
-0.0781 0.0562 -0.0503 0.0055 0.0858 0.0042 0.0018 -0.0083 0.0148 0.0613 -0.9826 0.1060
-0.8106 0.4936 -0.0154 0.3000 -0.0155 0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0938 0.0000
-0.8106 0.4936 -0.0154 0.3000 -0.0155 0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0938 0.0000
print(matrix_format(m))
-0.8106 0.4936 -0.0154
0.3000 -0.0155 0.0000
0.0000 -0.0000 0.0000
Отредактировано izekia (Ноя. 8, 2016 19:11:48)
Офлайн