#!/usr/bin/env python
# coding: utf-8
#
# 線形代数(Linear algebra) scipy版
#
#
# cc by Shigeto R. Nishitani, 2017-2018
#
# * file: /Users/bob/python/doing_math_with_python/linear_algebra_scipy.ipynb
#
# pythonでの線形代数のlibraryはいくつもあるように見える．
# 混乱しがちなんで，それをいくつか分類してみた．
# 間違ってたら教えて．
#
# * numpy.linalg
# * scipy.linalg
# * sympy
#
# scipyはnumpyのwrapperっぽい．[NumpyとScipy](https://www.eidos.ic.i.u-tokyo.ac.jp/~tau/lecture/computational_physics/slide/numpy.pdf)によると，
# 1. NumPy ⊂ SciPy ということのようだ
# 1. numpy で提供されている機能はそのまま, scipy でも提供されている
# 1. なので scipy だけで押し通しても良さそうだが, 世の中の説明は numpy が主流なので, それに合わせて, 基本は numpy, scipy だけで 提供されている機能は scipy を使う
#
# と明言されている．これが一番混乱がなさそう．
#
# sympyは代数計算向きなんで，表記がだいぶ違う．[別章](./linear_algebra_sympy.ipynb)として提供．
#
#
#
# # Table of Contents
#
# # 出力精度の制御
#
# 行列の出力はとても醜くなる場合がある．
# ほぼ0なのにとか，次元が増えてとかで．そのときこいつが役立つ．
# In[2]:
import numpy as np
from pprint import pprint
import scipy.linalg as linalg
a = np.array([[2,5], [4,1]])
l,P = np.linalg.eig(a)
pprint(l)
pprint(P)
np.set_printoptions(precision=3, suppress=True)
pprint(l)
pprint(P)
# # 行列，ベクトルの生成
# In[3]:
import numpy as np
import scipy as sp
list_a = [1,2]
list_b = [3,4]
np_a = np.array([list_a, list_b])
print(np_a)
sp_a = sp.array([list_a, list_b])
print(sp_a)
np_m = np.matrix([list_a, list_b])
print(np_m)
# In[4]:
list_a = [1,2]
np_a_1 = np.array(list_a)
print(np_a_1)
np_v = np.array([list_a])
print(np_v)
# np.vectorはない．
# # 行列の形状(shape)
# In[8]:
print(np_v.shape)
print(np_m.shape)
# ## ゼロ行列，単位行列
# In[10]:
np.zeros((3,3))
# In[8]:
np.diag([1,2,3])
# In[9]:
np.identity(3)
# ## 転置(transpose)
# In[7]:
print(np_a.transpose())
print(sp_a.transpose())
print(np_m.transpose())
# In[9]:
np_v.transpose()
# In[12]:
print(np_v.shape)
print(np_v.transpose().shape)
# # 行列，ベクトルの演算
# ## ドット積
# $$
# \left[
# \begin {array}{cc}
# 1&2\\
# 3&4
# \end {array} \right]
# \left[
# \begin {array}{c}
# 1\\
# 2
# \end {array} \right]
# $$
#
# を計算したいときは，
# In[4]:
np.dot(np_a,np.transpose(np_v))
# In[5]:
np_v.dot(np_a)
# これだと，
# $$
# \left[
# \begin {array}{cc}
# 1&2
# \end {array} \right]
# \left[
# \begin {array}{cc}
# 1&2\\
# 3&4
# \end {array} \right]
# $$
#
# を計算することになる．
# ベクトルの距離
# $$
# \left[
# \begin {array}{cc}
# 1&2
# \end {array} \right]
# \left[
# \begin {array}{cc}
# 1\\
# 2
# \end {array} \right]
# $$
# は
# In[13]:
np_v.dot(np_v)
# ではダメで，
# In[14]:
np_v.dot(np_v.transpose())
# In[15]:
np_v.transpose().dot(np_v)
# と順番を間違うと悲惨．単に長さだけなら$||v||$の代わりにnp.linalg.norm(v)も使える．
# In[18]:
np.linalg.norm(np_v)**2
# ただし，normなんでそのほかの指数も指定できる．
# In[19]:
np.linalg.norm(np_v,1)**2
# ## 外積，outer, cross
# In[30]:
v1 = np.array([1,1,3])
v2 = np.array([1,2,-1])
# In[31]:
np.outer(v1,v2)
# In[32]:
np.cross(v1,v2)
# ## スカラー三重積
# In[14]:
v3 = np.array([-1,2,1])
np.dot(np.cross(v1,v2),v3)
# # 行列要素のとりだし，追加
# In[33]:
a=np.array([[1,2,3],[4,5,6],[7,8,9]])
# i行j列の要素の取り出し
a[1,1]
# In[34]:
a[:2,1:4]
# In[35]:
col_0 = a[:,0]
print(col_0)
# In[36]:
row_1 = a[1,:]
print(row_1)
# In[37]:
b = np.array([1,2,3])
np.vstack((a,b))
# hstackはa,bが同じ次元でないとダメとのお達し．拡大係数行列を手軽に作れるわけではなさそう．column_stackを使うとできた．
# In[38]:
np.column_stack((a,b))
# # 掃き出し, LU分解
# In[43]:
import scipy.linalg
a = np.array([[1,-1,-1],[1,-1,1],[1,1,-1]])
P, L, U = scipy.linalg.lu(a)
print(P)
print(L)
print(U)
# In[45]:
b = np.array([1,-1,-1])
ab = np.column_stack((a,b))
P, L, U = scipy.linalg.lu(ab)
print(P)
print(L)
print(U)
# ## 階数
# In[42]:
print(a)
print(np.linalg.matrix_rank(a))
print(np.rank(a)) #deprecatedなんでやめろって，
# 関数名とかlibの区分けに統一性がないよね．．．
# # 逆行列
# In[9]:
a = np.array([[1,2],[3,4]])
scipy.linalg.inv(a)
# ## 行列式
# In[8]:
c = np.array([[1,1],[1,2]])
print(np.linalg.det(c))
a = np.array([[1,2],[3,4]])
scipy.linalg.det(a)
# ## 連立方程式の解
# 逆行列を用いて，連立方程式の解を求めるには次の通り．
# 例えば，連立方程式が次のような場合，
# $$
# x - y -z = 1 \\
# x - y +z = -1 \\
# x + y -z = -1
# $$
# In[10]:
a=np.array([[1,-1,-1],[1,-1,1],[1,1,-1]])
b = np.array([1,-1,-1])
scipy.linalg.inv(a).dot(b)
# ## 課題
# 次の連立方程式を解け．
# $$
# \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[26]:
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]))
# # 固有値，固有ベクトル
# In[11]:
a = np.array([[1,1],[0,2]])
l, P = scipy.linalg.eig(a)
print(l)
print(P)
# ## ベクトルの規格化
# In[12]:
# eigの固有ベクトルは規格化されている
p_norm = np.linalg.norm(P[:,1])
p_norm
# 規格化はこいつをつかって．．．
# ## 対角化
#
# Pに固有ベクトルが入っているので，対角化は$P^{-1} A P$で
# In[13]:
np.dot(np.dot(np.linalg.inv(P),a),P)
#np.dot(np.dot(np.linalg.inv(P),a),P)
# ## 課題
# 行列
# $$
# A\, = \, \left[
# \begin {array}{ccc}
# 2&0&1\\
# 0&3&0\\
# 1&0&2
# \end {array} \right]
# $$
#
# を対角化する変換行列$P$を求め，対角化せよ．
# In[30]:
import numpy as np
A = np.array([[2,0,1],[0,3,0],[1,0,2]])
l, P = np.linalg.eig(A)
print(l)
print(P)
np.dot(np.dot(np.linalg.inv(P),A),P)
# In[31]:
import numpy as np
from math import pi
np.linspace(0,2*pi,9)
# In[ ]: