#!/usr/bin/env python
# coding: utf-8
#
# 線形代数(Linear algebra)
#
#
# cc by Shigeto R. Nishitani, 2017-10-31
#
# * file: /Users/bob/python/doing_math_with_python/linear_algebra.ipynb
# Table of Contents
#
# # sympy
#
# * [sympy tutorial](http://docs.sympy.org/latest/tutorial/index.html)
# * [module references](http://docs.sympy.org/latest/modules/matrices/matrices.html) LUdecompositionの記述あり．
# # 行列，ベクトルの生成
# ## 直接的に
# In[1]:
from sympy import *
#init_session()
#init_printing(use_unicode=True)
init_printing()
list_a = [1,2]
list_b = [3,4]
A=Matrix([list_a, list_b])
A
# In[2]:
x,y = symbols('x,y')
v1=Matrix([x,y])
v1
# ## ゼロ行列，単位行列
# In[3]:
zeros(2,3)
# In[4]:
ones(2,3)
# In[5]:
eye(3)
# In[6]:
diag(1,2,3)
# ## 行列要素のとりだし，追加
# In[7]:
A=Matrix([[1,2,3],[4,5,6],[7,8,9]])
# i行j列の要素の取り出し
A[1,1]
# In[8]:
A[:2,1:4]
# In[9]:
col_0 = A[:,0]
col_0
# In[10]:
row_1 = A[1,:]
row_1
# numpyではcolumn_stackを使って拡大係数行列を作った．
# sympyではrow_insert, col_insertを使う．
# In[11]:
b = Matrix([1,2,3])
A.row_insert(3,b.T)
# In[12]:
A.col_insert(3,b)
# # 行列・ベクトル間の単純な演算
# ## 次元の確認
# In[13]:
A.shape
# In[14]:
print(v1.shape)
print(v1.T.shape)
# ## 転置(transpose)
# In[15]:
A.T
# In[16]:
v1.T
# ## ドット積(内積)
# numpyと違ってdot関数を使うと怒られる．
#
# "\*" を使います．
# In[17]:
A = Matrix([[1,2],[3,4]])
A*v1
# In[18]:
v1.T*v1 #v1*v1ではerror
# In[19]:
v1.T*A
# In[20]:
v1.T*v1
# In[21]:
v1*v1.T
# ## 外積，outer, cross
#
# outer(v1,v2)がSympy1.1.1では見つからない．
# crossは使える．
# In[22]:
v1 = Matrix([1,1,3])
v2 = Matrix([1,2,-1])
# In[23]:
v1.outer(v2)
# In[24]:
v1.cross(v2)
# ## スカラー三重積
# In[25]:
v3 = Matrix([-1,2,1])
v3.dot(v1.cross(v2))
# # 行列操作
# ## 掃き出し, LU分解
# In[26]:
A=Matrix([[1,2,3],[4,5,6],[5,7,9]])
print(A.nullspace())
Aex = A.col_insert(3,b)
pprint(Aex)
Aex.rref()
# rrefはreduced row echelon formを取り出してくれる．これはGaussの消去法で下んとこまで行って，さらに上んとこへ戻ってきたやつ．拡張行列もそのままやってくれるはず．もうすこし，試す必要があるが，一応MapleのLUDecompositionに近い結果を与えてくれる．
# In[27]:
A=Matrix([[1,2,3],[4,5,6],[7,9,8]])
L, U, P = A.LUdecomposition()
L, U, P
# ## 階数, rank
# In[28]:
print(A)
print(A.rank())
# ## 行列式, determinant
# In[29]:
c = Matrix([[1,1],[1,2]])
print(c.det())
a = Matrix([[1,2],[3,4]])
a.det()
# ### 課題
# 次の行列式の計算を確認せよ．
# $$
# \left|
# \begin {array}{cccc}
# 1+a & 1 & 1 & 1\\
# 1 & 1+b & 1 & 1\\
# 1 & 1 & 1+c & 1\\
# 1 & 1 & 1 & 1+d
# \end {array} \right|
# = abcd \left(1+\sum_{i=a}^d \frac{1}{i}\right)
# $$
#
# In[30]:
from sympy import *
a,b,c,d,x,y,z = symbols('a b c d x y z')
A = Matrix([[1+a,1,1,1],[1,1+b,1,1],[1,1,1+c,1],[1,1,1,1+d]]
)
pprint(A)
simplify(A.det())
# ## 逆行列
# In[31]:
A = Matrix([[1,2],[3,4]])
A.inv()
# ### 課題
# 次の連立方程式を解け．
# $$
# \left\{
# \begin {array}{cl}
# x+y+z&=0\\
# ax+by+cz&=0\\
# bcx+cay+abz&=(a-b)(b-c)(c-a)
# \end {array} \right.
# $$
# In[32]:
from pprint import pprint
from sympy import *
a,b,c,x,y,z = symbols('a b c x y z')
A = Matrix([[1,1,1],[a,b,c],[b*c,c*a,a*b]])
bb = Matrix([0,0,(a-b)*(b-c)*(c-a)])
print(A)
print(bb)
Ainv = A.inv()
res = Ainv * bb
pprint(simplify(res[0]))
pprint(simplify(res[1]))
pprint(simplify(res[2]))
# # 固有値，固有ベクトル
#
# sympyでは固有値・固有ベクトルは，eigenvects()で求められます．
# 戻り値は，
# > [(固有値, 重複度, [固有ベクトル]), ...]
#
# です．
# In[33]:
A = Matrix([[1,2],[2,1]])
A.eigenvects()
# ## ベクトルの規格化
#
# eigenvectsの固有ベクトルは規格化されていません．norm()を使うとかんたんに規格化してくれます．
# In[34]:
v1=Matrix(A.eigenvects()[0][2])
v2=Matrix(A.eigenvects()[1][2])
simplify(v1.T*v1)
# In[35]:
v_norm = v1.norm()
v_norm
# In[36]:
v1n = v1/v_norm
v1n
# ## 対角化
#
# ちょいとわかりやすい行列を作成してみてみます．
# In[37]:
D = diag(1,2)
P = Matrix([[1,0],[1,1]]).T
A = P*D*P.inv()
A, P, D
# In[38]:
P,D = A.diagonalize()
P.inv()*A*P
# In[39]:
eigs = A.eigenvects()
v1 = eigs[0][-1][0]
v2 = eigs[1][-1][0]
Matrix([v1,v2]).reshape(2,2).T
# 対角化行列$P$をeigenvectsから作ろうとして苦労した．まず，固有ベクトルの取り出しが自動的にできない．重複しているときには，固有ベクトルは二つ出てくる．さらに，タプルからの取り出しがこれでいいいのか．．．さらに，ベクトルから行列を作るときに横置きのができない．うううんん．他のやり方をしたでやってみる．上では仕方がないので，一旦4x1のMatrixにしてそのあとreshape，さらにTransposeで．．．
# In[40]:
v1.col_insert(1,v2)
# numpyにはravel()というフラット化の関数があるがsympyにはない．
#
# 結論：3つほど上でやっているdiagonalizeで取り出しましょう．重複があっても取り出してくれるんで．
# ### 課題
# 行列
# $$
# A\, := \, \left[
# \begin {array}{ccc}
# 2&0&1\\
# 0&3&0\\
# 1&0&2
# \end {array} \right]
# $$
#
# を対角化する変換行列$P$を求め，対角化せよ．
# In[41]:
A = Matrix([[2,0,1],[0,3,0],[1,0,2]])
evs = A.eigenvects()
evs
#
# In[42]:
P,D = A.diagonalize()
P,D
# In[43]:
P.inv()*A*P
# In[ ]:
# In[ ]: